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Abstract 

We present not a literature review but a description, as detailed and consis¬ 
tent as possible, of two analytic models of disk accretion onto a rotating black 
hole: a standard relativistic disk and a twisted relativistic disk. Although one of 
these models is much older than the other, both are of topical current interest 
for black hole studies. The way the exposition is presented, the reader with 
only a limited knowledge of general relativity and relativistic hydrodynamics 
can — with little or no use of additional sources - gain good insight into many 
technical details lacking in the original papers. 
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Chapter 1 

Relativistic standard 
accretion disk 


In Chapter 1, the model of a standard accretion disk around a rotating black 
hole is presented with general relativity effects rigorously taken into account. 
This model was first described in paper [T] and since then has been used in 
many studies to obtain a convincing evidence of the existence of black holes 
in both stellar binary systems and active galactic nuclei. It remains topical 
because a full account for general relativistic properties of matter motion in the 
disk and generation of its emission allows the position of the inner disk radius 
and hence the black hole spin to be inferred from observations. Observational 
appearances of relativistic disks were modeled for the first time in [5] (see [3], 
[3] for comprehensive reviews, references therein and their citations list). In 
addition, the standard accretion disk underlies more complicated theories of 
warped (twisted) accretion disks which are formed when the accreting matter 
moves outside the equatorial plane of a rotating black hole. Such a theory is 
presented in [S] and is discussed in more detail in Chapter 2. 

Everywhere below, the natural units G = c = 1 are used. If the mass is 
measured in units of the black hole mass, M, the unit of length is half of the 
Schwarzschild gravitational radius, Rg/‘2., such that Rg/2 = GM/(? = 1 and 
the unit of time is the light crossing time of the unit of length. 

In addition, Latin indices i,j,k... taking values from 0 to 3 are used to 
denote components of vectors, with the zero component standing for the time 
coordinate. Also, wherever needed the Einstein summation convention is used. 


1.1 Space-time near rotating black holes 

1.1.1 Kerr metric 

Properties of spacetime near a rotating black hole are described by axially sym¬ 
metric and stationary metric of the form (see, e.g., [B], paragraph 4.2): 

— ujdt)^ — ^dR? — g^dO^, (I.l) 
A 
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where the signature (1,-1,—1,-1) is chosen and the coefficients are 

= {R^ + a^ f - a^Asin^e, 

g'^ = R^ cos^ 9 and A = R^ — 2R -\- a? 

The coordinates {t, </>, i?, 9} are called the Boyer-Lindquist coordinates. Far 
away from the gravitating body, the spatial part of these coordinates in the 
limit of the zero black hole spin parameter, a, transits into the usual spherical 
coordinates, where (j) is the azimuthal angle. For non-zero a, it transits into gen¬ 
eralized spherical coordinates in which the surfaces of constant radial distance, 
R = const, represent spheroids with the aspect ratio R/{R? + 

The space-time described by (HH), is axially symmetric with respect to the 
line 9 — Q, called the black hole rotation axis; the plane corresponding to 9 = 
7r/2, is called the black hole equatorial axis. 

In (11.11) an important quantity appears: 


UJ = 


2aR 

'W’ 


( 1 . 2 ) 


which has the dimension of frequency. This is the angular velocity that every 
freely moving observer, without exception, acquires in the direction of the black 
hole rotation. 

As is described in the literature on the structure of rotating black holes (see, 
e-g-, [Z], paragraph 58), metric (11.111 has several special hypersurfaces, including 
the event horizon and the ergosphere. However, as we show below, for the 
astrophysical problem under consideration, of most importance is the dynamics 
of free circular motion of particles in the equatorial plane of the gravitating 
body. This motion has additional features in comparison to the corresponding 
Newtonian problem. Note also, that in Chapter 2 weakly elliptical orbits slightly 
inclined to the equatorial plane are considered. 

We consider a standard, and hence geometrically thin, accretion disk. Such 
a disk is basically flat. By definition, this is a stationary flow of matter with 
mirror symmetry with respect to its middle plane and axial symmetry with 
respect to the line perpendicular to this plane. Clearly, such a model flow can 
be described by dynamical equations in axially symmetric metric only if the 
disk symmetry plane coincides with the equatorial plane of the black hole. To 
tackle the problem, only the form of metric near the plane 9 = 12 is sufficient. 
Passing to cylindrical coordinates using the standard transformation 


r = i?sin0, z = Rcos9, 


all metric coefficients gik in (HU can be expanded in a power series in small 
ratio zjr ^ 1. For geometrically thin disks, the corrections to gtk due to non- 
equatorial motion up to (z.jrf' are sufficient. Indeed, one of the basic equations 
describing the disk, namely, the projection of the relativistic analog of the Euler 
equation onto the direction normal to the disk plane, must be odd with respect 
to the coordinate reflection z —>■ —z due to the mirror symmetry of the disk. 
This means that in its series expansion in (z/r) only odd powers of (z/r) must 
be present. By the main assumption on the smallness of (z:/r) only the first term 
in this expansion should be kept. This, in turn, corresponds to series expansion 
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of gik up to quadratic terms, since only first derivatives of gik, characterizing 
the ’strength’ of the gravitational field, enter dynamical equations. 

Note, however, that hydrodynamic equations also contain a second covariant 
derivative of the velocity field (see below), and hence the final expressions can 
involve second derivatives of gik with respect to z, which may seem to require 
that we keep terms of the order of {z/r)^ in gik- But this is not required, 
because, as follows from the explicit form of the stress-energy tensor, such terms 
can appear only when multiplied by some of the viscous coefficients, which in 
turn cannot be greater than of the order of (z/r) being proportional to the 
characteristic mixing length in the fluid. The latter is initially assumed to be 
less than the disk thickness. 

As regards other equations, namely (see below): two projections of the rel¬ 
ativistic analog of the Euler equation on to the disk plane, the energy balance 
equation and the rest energy conservation law - the same symmetry considera¬ 
tions imply that they are even under the coordinate reflection z —>■ —z; therefore, 
the leading term is of the zeroth order in (z/r) in the metric expansion. 

Using these expansions and expressions for the coordinate differentials. 


dR = 



dr H—dz, 
r 


de = 


z dr 



dz 

5 

r 


we find the metric in the following form (see also 0 ): 


ds^ = 


1 - 


I 1-b 

2a 
r 


2a^ 


dv— 

2-il3 


2 I 2 

r -I- a 


2a^ 




r 

dtd(j)— 


l + --b 


r^D 


--4-^13-^ 


6 2a^ 

TT 


— ^ 1 drdz— 


1-b 


dr'^ 


D rD \ r 


2a^ 


r^D 


dz^, ( 1 . 3 ) 


where the notation 


£1 = 1- 


is introduced. Below, we also use (with a few exceptions) the notations in¬ 
troduced in the original paper by Novikov and Thorne [I] for the relativistic 
correction coefficients. 

Finally, the inverse of the matrix g'^^ corresponding to double-contravariant 
tensor has the form: 


5*'= = 


{gttg^^ - gtcf,) x 


-9t^ 


9t(t> 

-9tt 


{grrgzz 9rz) X 


gzz 

-9rz 


-grz 

~9rr 

( 1 . 4 ) 


d(f^+ 
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1.1.2 Circular equatorial geodesics 

The expression for circular equatorial geodesics can be conveniently found from 
the extremum condition for the distance along them. Here we follow the exposi¬ 
tion from [9], (see paragraphs 13.10 and 13.13). Indeed, for time-like trajectories 
the functional 


5 = Lds = 


dx^ dx^ , 
Qik q ds, 
ds ds 


should be minimal, which is equivalent to the Euler-Lagrange equations for L\ 


A. (—'] 

ds \dx^) 



(1.5) 


where Ug = dx'^/ds = xd is the four-velocity in the Boyer-Lindquist coordi¬ 
nates. As L does not explicitly depend on t and (j), the following quantities are 
conserved: 

9uUl = k, 

where k and h have the meaning of the time and azimuthal covariant velocity 
components, respectively. 

In explicit form, using the components gik from (irnD at z = 0, we find: 



2a • , 

t -I-(/) = k, 

r 


( 1 . 6 ) 



r 



= —h. 


(1.7) 


We temporarily assume that the motion is not necessarily circular and 
Ug ^ 0. Instead of the r-component of the Euler-Lagrange equations, it is 
more convenient to use the condition of the normalization of the four-velocity 
of particles with non-zero mass: 


- 2g*'^kh + g'^'^h'^ + = 1 . 

This yields the following equation for k and h: 

P — 1 

Y + r.,/M = ^-. 

where we introduce the effective potential 

1 h'^ — a^{k^ — l) {h — akY 


( 1 . 8 ) 


(1.9) 


( 1 . 10 ) 


The conditions for circular motion include, first, r = 0 and, second, r = 0 
(for the particle to stay in a circular orbit). The latter condition is equivalent 
to the vanishing of the derivative of 14// with respect to r: 


1-b —^- - - 

r 


3{h — akY 


= 0 . 


( 1 . 11 ) 
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Equation (HU with f = 0 and equation (11.111) allow us to determine k and 
h as functions of r and then, using (11.61) and (11.71) . to find U* and U^. 

To solve the first problem, let us introduce the new variable fj, = h — ak 
and, to facilitate manipulations, make the change u = 1/r. Then equation (11.911 
taken with r = 0 and equation ((1.111) yield the following equation for /i: 

m^[(3m— 1)^ — — 2u[{3u — l)(o^u ~ 1) ~ 2ua^{u— 1)]^*^ + (art — 1)^ = 0. 

( 1 . 12 ) 

The solution of (11.12(1 for a stable circular prograde orbit has the form: 


M = 


ly/u — 1 


[it(l — 3it + 2au^/^)]i/2 


(1.13) 


Using (11.1311 and (11.91) taken at r = 0, we find the constants h and k, as well 
as the components U*: 


Ug* = B, Ug'>’ = (r3C')-^/^ t/g’' = Ug^ = 0, (1.14) 


where 




(1.15) 


It is easy to check that modulus of this vector is equal to unity: 


g^k Ug^Ug^ = 1 


The angular velocity measured by the clock of an infinite observer (who 
measures the coordinate time t ), corresponding to such a motion, is 

U = ^ = r-3/2i?-i (1.16) 

dt ^ ’ 

It follows that in the Schwarzschild case this value exactly coincides with the 
Keplerian angular velocity. 


1.1.3 Radius of the innermost (marginally) stable orbit 


This is determined by the condition that the stable circular motion is no longer 
possible when the minimum of the function Veff{r,h{rc),k(rc)) disappears at 
r = Tc, where Tc is the radius of a circular orbit. This is equivalent to the 
condition 


d^Veff 

dr'^ 


r—Tc 


= 0 , 


which leads to the quartic equation 


— 6z^ + 8az — 3a^ = 0, 


(1.17) 


where z = r^/^. 

Using the Ferrari method (see, e.g., [10(1. we write the corresponding auxil¬ 
iary cubic equation: 


y3 - 12?/2 + 12(3-t a2)y - 64a^ = 0. (1.18) 
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The real root of equation (11.181) is related to the Cardano solution to the corre¬ 
sponding incomplete cubic equation and is given by 

2/1 = -2(1 - a 2 )i/ 3 [(i + a)i /3 + (i _ a)i/ 3 ] + 4, (419) 

Next, having j/i, it is possible to use the Ferrari solution to write the 
quadratic equation that gives two real roots of (11.171) : 

/ + VinP + \ (^-6 + 2/1 - = 0 ( 1 . 20 ) 

The larger root of (11.201) . pi, determines the boundary of the stable circular 
motion of test particles in the equatorial plane, which we denote as r = r^s- 
Thus, 

rms=pl=3 + - (-J/?/4 -b (1.21) 

\/yi 

It is easy to check that the result (11.211) coincides with the expression presented 
in [a, (see formula (15k) therein), taking into account that the auxiliary values 
Zi .2 in [llj take the form Zi = 3 — yi/2 and Z 2 = ^aj^Jy\ for a > 0 in our 
notation. 

In the case of the Schwarzschild metric, o = 0, we recover the well-known 
result that the circular motion becomes unstable for r < 6, i.e. at distances 
smaller than three gravitational radii from the black hole. For slow rotation, 
1 ^ a > 0, we have r^s ~ 6 — 4-\/6a/3, and hence the zone of stable motion 
shifts closer towards the event horizon. In the limit case a = 1 we find r^s = !> 
i.e. the marginally stable circular orbit coincides with the gravitational radius 
of the extreme-spin black hole. 

During accretion, gas elements in the disk slowly approach r^s by loosing 
their angular momentum due to the action of viscous forces. Once the gas 
elements fall into the region with r < rms, due to instability of the circular 
motion they need not lose the angular momentum any more to approach the 
black hole. This means that the matter falls freely inside r^s, and the standard 
accretion disk model assumes that r^s is the inner disk radius. 


1.2 Choice of the reference frame 

1.2.1 Bases in general relativity 

Mechanical laws, formulated in the form of vector equations, can be written in 
the symbolic form irrespective to any observer or any reference frame. But to 
represent some physical quantity describing a natural phenomenon in the form 
of a set of numerical values, the measurement procedure should be specified. 
In Newtonian mechanics, this means that the observer introduces a coordinate 
system, and then at each point of space he/she arbitrarily constructs three basis 
vectors. The coordinate system and the basis vectors can evolve in time, which is 
the same at all points. A tool measuring time, together with the coordinate grid 
and a vector basis, form a reference frame, in which any physical value (scalar, 
vector or tensor) can be measured, i.e. it can be represented by a number or a 
set of numbers. 
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The situation in relativistic mechanics is different: because it is not possible 
to consider the time independently, it becomes the fourth component of the 
space-time continuum. Therefore, the choice of the reference frame reduces to 
the construction of a coordinate system and four basis vectors determined at 
each point of space-time. In general, this set of basis vectors is usually referred 
to as a tetrad. There is no a universal observer any more; instead, a set of 
observers moving along a certain family of world lines is considered. If one of 
the tetrad orts, conventionally corresponding to the time direction, is tangent 
at each point to these world lines, the tetrad is said to be ’transported’ by the 
observers. The last statement can be easily understood, because in such a basis, 
the four-velocity of each observer at any time has a non-zero projection only 
on the ’time’ ort, in other words, the observers are at rest in this basis, i.e. 
transport it with them. 

Coordinate representation 

Thus, the choice of the coordinate system and the choice of the tetrad are 
independent procedures. Nevertheless, if there is a coordinate system, x®, the 
tetrad is frequently chosen in such a way that each basis vector, e^, is tangent 
to the corresponding coordinate line. Here the moduli of orts of these so-called 
coordinate bases are chosen such that their pairwise scalar products are equal 
to the corresponding metric coefficients: 

We recall that in differential geometry (see paragraphs 3.1-3.4 i), such 
coordinate orts are introduced as objects isomorphic to the partial derivatives 
of an arbitrary scalar function on the manifold with respect to the coordinates. 



Any (tangent) vector is a linear combination of the coordinate orts, and the 
components of this linear combination are called contravariant components of 
the vector. 

In addition to ej the so-called dual basis, e‘, is introduced with the orts 
defined as 

(ei-ej) = d,^ (1.24) 

where <5/ is the Kronecker symbol. The condition (11.241) implies that each ort 
of the dual basis has unit projection on the corresponding ort of the coordinate 
basis and is orthogonal to all other orts of the coordinate basis. 

The dual coordinate orts, in turn, are introduced as objects isomorphic to 
the coordinate differentials, 

= dx\ (1.25) 

Next, if we use the fact that any tangent vector A can be alternatively 
presented as a linear combination of dual coordinate orts, whose coefficients are 
referred to as covariant vector components, we obtain the well-known rule of 
lowering vector indices: 

Ak = Aj(e* • Gfe) = (A^e* • efc) = (A'e^ • e^) = A\e, ■ Sk) = A’-gik, (1.26) 
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In a similar way it is easy to show that if we introduce the notation = (e®-e^), 
then due to the duality of bases, the matrix is inverse to the matrix gik, and 
the rule of raising of vector indices holds. Similar representation in coordinate 
bases can be extended to a more general case of tensors. 

Tetrad representation 

In this and subsequent sections we mostly follow the exposition given in para¬ 
graph 7 of [12]. Assume that now we want to project the same vectors and 
tensors on an arbitrary tetrad defined by relations 

e(o) = e(o) (1-27) 

where e(a)® are some functions of coordinates, and the indices labeling the tetrad 
orts are in parentheses. 

From the duality condition (|1.24p we can introduce the dual tetrad: 

= e(“),dx®, (1.28) 

where is the matrix inverse to e(a)b 

In these matrices, there are two kinds of indices: coordinate and tetrad ones. 
The coordinate indices can be lowered or raised using the metric (HI). We can 
impose an additional constraint on the tetrad: 

e(a)*e(6)^ = ?7(a)(b), (1.29) 

where 

(1.30) 

are mutually inverse matrices and r]{c){b) = diag{l, —1, —1, —1) is the Minkowski 
matrix. In other words, we require that the original and dual tetrads be or¬ 
thonormal in four-dimensional pseudo-Euclidean space. 

Using above relations, it is straightforward to show that 

( 1 . 31 ) 

and therefore the following alternative expression for the interval squared holds: 

ds^ = dx^) dx’") = 'n(a)ib) (1.32) 

which is useful below. 

Note that the values in parentheses in the right-hand side of (11.321) can be 
considered as infinitesimal shifts along the corresponding orts of the tetrad; 
therefore, in the introduced tetrad representation with an orthonormal tetrad, 
the square of interval takes exactly the form as in the Minkowski space-time of 
special relativity. Similarly, expressions (11.271) can be thought of as directional 
derivatives along the tetrad’s orts, and these have exactly the form that the 
usual partial derivatives with respect to coordinates in the coordinate basis 
take when changing from the coordinate basis to the tetrad one. 

Using the definitions and relations given above, it is easy to see how the 
tetrad components of vectors are expressed through the coordinate components. 
Tetrad components of a vector are written as 

A(,) = e®„)A„ = e(“)A® = (1.33) 
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Conversely, 


A . — 4 , , Ai — pi 4(“) 

Similar expressions can be written for tensor of any valence. For example, 
for a two-covariant tensor, we have 

T{a}{b) = ~ ^{a)"^i{b), 

and conversely, 

(1.34) 

Note in conclusion that relations (jl.ddp and isomorphism (11.251) can be used 
to find contravariant components of the four-velocity in the tetrad representa¬ 
tion: 


= —. (1.35) 

as 

That is, it is again a unit tangent vector along the world line, but now its com¬ 
ponents are given by small shifts along the corresponding orts of the dual basis. 
Using (11.3311 it is easy to find the relation between the conventional coordinate 
components of four-velocity, U® = dx'^/ds, and its tetrad components. Covari¬ 
ant tetrad components are derived from contravariant ones using the standard 
rule in special relativity: lowering a spatial index is equivalent to changing the 
sign of the corresponding component. 


Covariant derivative in tetrad representation 

Let us calculate the directional derivative along a tetrad ort from a contravariant 
component of a vector: 


A 


a 


(a),(b) dx 


a 


'(a) 


-(b) 


■Ake\ 


a);2J ’ 


(1.36) 


where the semicolon denotes the usual covariant derivative in the coordinate 
basis. 

Expression (11.361) can be recast to the form: 



4l-(a),(h) — + ^{a)k-i^{b)^(c)-^^ 

(1.37) 

whence 

^\a)-^rd^(b) = ^(a),(b) ~ 7(c)(a)(b)^^ K 

(1.38) 

where 7(c) (a) (6) 

are the so-called Ricci rotation coefficients. 



7(a)(b)(c) = e(b)fc;*e(^)e^„) 

(1.39) 


An important point is that for orthonormal bases satisfying (ll.29|l . the coeffi¬ 
cients 7(a)(b)(c) are antisymmetric in first two indices. Indeed, 


0 — {V(b)(a)),i — (e(b)fce^a));i — ^(b)k-ie\a) + ^(b)ke\a)-i — ^(b)k-ie'la) + 

Comparing this relation with (11.391) proves the stated property of the Ricci 
coefficients. 
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Finally, let us discuss one more useful property of coefficients (11.391) : to cal¬ 
culate these coefficients only partial derivatives of the components of the tetrad 
basis orts are needed, and therefore the Christoffel symbols are not required. 
Indeed, we consider auxiliary combinations 

^(a)(f))(c) = ^(b)i,jl^(a)^lc) ~ ®(a)®(c)]> (1-40) 

and rewrite them in the form 

^{a)(b)(c) = ®(a)®(c) ~ ^(b)j,i]’ (1-41) 

In the last expression, the ordinary partial derivatives can be substituted by co¬ 
variant ones, because the additional terms with Christoffel symbols are symmet¬ 
ric in i,j. Then expression (11.4111 is equal to the difference 7(a)(6)(c) ~ 7(c)(&)(a)- 
But in such a case, 

7(a)(b)(c) = l/2[A(a)(b)(c) + A(c)(a)(b) “ A(b)(c)(a)] (1-42) 

and, using (I1.40L it is possible to calculate the Ricci rotation coefficients by 
taking partial derivatives of components of the tetrad basis orts. 

We now consider formula (11.381) . The left-hand side represents simply the 
projection on tetrad basis of a rank-2 covariant tensor obtained by taking the 
derivative of some vector field. Therefore, this combination has meaning of the 
covariant derivative of a vector taken in a non-coordinate basis. 

Next, the right-hand side of (11.381) has exactly the same form as the covari¬ 
ant derivative in a coordinate basis, with the only difference that it involves 
tetrad indices (which can be raised or lowered, including for 7(o)(b)(c); using 
the Minkowski metric). It can be shown that the same holds for contravariant 
components of a vector field and for tensor fields in general. 

Thus, as the components of a covariant derivative in a tetrad basis has the 
same form as in a coordinate basis, it is convenient to use the same notations 
and terms that are used in the coordinate basis. In particular, the Ricci rotation 
coefficients are simply referred to as the connection symbols in a given basis. 
We emphasize once again that they should not be confused with the Christof¬ 
fel symbols, which represent another limit case of connection coefficients in a 
coordinate basis and have different index symmetry. 

I. 2.2 Tetrad transported by rotating observers 

We construct a tetrad basis related at each point of space-time to observers 
moving around a black hole in equatorial circular orbits with angular velocity 

n. 

Strongly at z = 0 this is the free motion along geodesics found in Section 

II. 1.21 However, for a small deviation from the equatorial plane, such a motion, 
corresponding to the constant z, is possible only if there is some external sup¬ 
porting force; in the case of a gas disk, for example, this force is due to the 
pressure gradient. 

To start the construction, we direct the time ort of the tetrad along the 
world line under discussion. Using the four-vector of the geodesic found in 
Section fl. 1.21 write it in the form 

e(,) = (t7/ + Zo)| + t// 
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where we add the correction factor Zfj^zjr) to the time coordinate component, 
since the modulus of the vector e(() should be equal to unity away from the 
equatorial plane as well, whereas the vector Ug itself is unitary only at z = 0. 
Clearly, with account for this correction, e(t) would correspond to the four- 
velocity of the real motion. Calculation of the modulus of the vector 0 ( 4 ) in 
metric USD shows that it is equal to unity under the following condition: 


7 rV ^ 

° Vr/ 2rGCi/2’ 

where we introduce the relativistic correction coefficients 



(1.43) 


4a 3a^ 
3/2 ^2 


(1.44) 


Thus, the ort e^t) is transported by the observer rotating around the black 
hole with a frequency equal to the (/)-component of ep), which is independent 
of z. This frequency corresponds to the free circular motion in the equatorial 
plane of the black hole, and rotation occurs in planes of constant z. 

We now calculate the time ort of the dual basis. According to the convention 
rule for raising and lowering coordinate indices, we have 


— (Ug^gtt + Ug'^gt^) dt + {Ug^g^t + Ug^^g^^) d(f>. 


Next, we consider the part of metric (ESI) containing differentials dr and dz. 
It can be rewritten in the form (see the result (11.321) in the previous section): 


dSrz 




where 

eW = ^ 

/ a 2\l/2 

are the radial and vertical orts of the dual basis, respectively. The coordinate 
components of vectors and satisfy the orthonormality condition 

(11.291) . as can be easily checked by direct substitution. 

The orthonormality condition for a tetrad can now be used to determine the 
fourth ort corresponding to the azimuthal direction. 

From three orthonormality conditions for three already known vectors, we 
obtain that for these conditions to be consistent, the following relation should 
hold: 


eW, = = 0 , 

and the time and azimuthal components should be related as 


pW, --pW, .fiE 

C 0 — c t , 
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Finally, the normalization condition for yields a quadratic equation for 
and the sign of the solution is dictated by the additional requirement of 
the choice of a right-hand triple of space orts of the tetrad. 

We thus obtain the dual tetrad basis with the leading corrections in (z/r) 
due to out-of-equatorial-plane motion in the form 


(*) — 




C-i/2|r./2F+(i)"^Z,| #, (L45) 


e(« = _ 


rC J 
1/2 


1/2 


1 f z\'^ 1 — a/r 


(;) 


2 Vr 


r3/2 


dt + 


„ fDy^^ 1 /zy , B H [D\ 

+51;) ^ g i c j 


1 / 2 ' 


dcj), (1.46) 




(-) = ( 1 + 4 ) dz. 


(1.48) 


To obtain the original basis, which we use to write equations of motion, it 
suffices to calculate the inverse to the matrix which yields 


"(t) 


= 


B- 


0\2 H 


(;) 


rJ 2rG 


| + (ec)-'G|, (1.49) 


6 ( 0 ) = 


(rC'£>)i /2 


O 


G 


d 


dt \r(£»C')i/2 


+ O 


d<j) 




1 (zy Z2 \ _d_ 

9r ’ 


2 9 , . 9 z^\ d 


®(-) = -12(2-0 A) + 1 - 


, (1.50) 

(1.51) 

(1.52) 


The following notations for relativistic correction coefficients are introduced 
in the expressions for the original and dual bases: 


2a 

j,3/2 j,2 


(1.53) 


5 a 3a 3a^ 2a^ 

^ — j,l/2 j,3/2 ^7/2 


(1.54) 


13 



(1.55) 



Here, we omit terms ~ 0{z‘^jr^) in the expression for the azimuthal ort of 
the original basis due to their complexity; in addition, as we will see below, 
these terms are not required in the standard accretion disk model. 

For the reader’s convenience, we here preserve the notations introduced in 
paper [T] for the coefficients B, C, D, F, G, but use the standard style of Latin 
letters, which is more familiar to the reader. In addition, the coefficient H is 
equivalent to the coefficient C introduced in [8]. We also draw attention to the 
fact that two other coefficients, A and B, introduced in the same paper, are 
equivalent to our coefficients D and C, respectively. It can be checked that 
the original and dual bases presented in [T] coincide with bases derived here at 
z = 0. 

Using formulas (11.45111.48)) and ()1.35ll . it is easy to deduce that solution 
(11.141) indeed yields 1/^“^ = (1,0, 0,0) in the equatorial plane. 

The connection coefficients 

Using ()1.40l) and then (11.421) and knowing the matrices of the original and dual 
bases given above, we can calculate the connection coefficients 7(a)(f))(c)- 

Of the 64 coefficients, 16 are equal to zero due to anti-symmetry of 7(a)(f))(c) 
in the first two indices. For the same reason, of the other coefficients, only half 
(i.e. 24) have to be found. Because we are interested in the region near the 
equatorial plane of the black hole, it makes sense to separate these coefficients 
into two groups: those that are ^ {z/r)^ in the leading order, and those that 
are proportional to the first power of {zjr). As mentioned in Section fl. 1.11 the 
latter coefficients must appear in the vertical projection of the relativistic Euler 
equation, while the former emerge in other equations. 

It can be shown that 

1) if there is no index (z) among the indices of 7(o)(b)(c); then 7(a)(f))(c) ^ {z./rY', 

2) if only one such index is present, then 7(o)(6)(c) ^ (z/r), and, finally, 

3) if two indices (z) appear in j(a){b){c), then this coefficient is of the second 
order in [zjr). 

Indeed, we examine formula (11.401) . Here the brackets contain the original 
basis components, which are summed with the coordinate derivatives of the dual 
basis components (the raising of a tetrad index can only change the sign of the 
component). 

In case (I) (a),(&),(c) ^ {z). As the (t) —, (</>)— and (r)—orts of the original 
basis have no z-component, only terms which do not contain derivatives with 
respect to z of the dual basis components and have no z-component of the dual 
(r)—ort make a non-zero contribution to 7(a)(b)(c)- Only in these two cases can 
the contribution ^ (z/r) appear, and hence we prove the statement (1). 

Now, in (11.401) let (6) = (z) and (a), (c) (z). Then the non-zero contri¬ 

bution can only be due to terms containing the z-component of the (t)—, (</>) — 
and (r)—orts of the original basis which are absent. Therefore, to check case (2) 
we should consider only the variant when in (11.401) (a) = (z) or (c) = (z). Here, 
the terms containing separately either r— or z—components of the (z)—ort of 
the original basis contribute. In the first variant, the proportionality to ^ (z/r) 
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is due to exactly the component while in the second, it is due to the deriva¬ 
tive with respect to 0 of one of the dual basis components that is always even 
function of z, as can be easily verified. 

We leave it to the reader to prove the statement (3). 

The counting shows that there must be 9 connection coefficients without the 
index (z), and hence an even function of z, and 12 coefficients with the index (z) 
and hence odd function of z. The calculation indicates that only 4 coefficients 
of the first type are non-zero, namely: 

1 H _3,2 

7(t)(0)(r) = 7(i)(r)(0) =-r- (1.56) 

'ywir)(t) = —r ^ l{<j>)(r){4>) = ^ ) (1-57) 

To compute coefficients (11.5611.571) it suffices to use bases taken without 
corrections in z. When constructing the standard disk model, the following 
facts are also important. First, the direct calculation shows that another 5 
connection coefficients of this type are zero through corrections of the order 
of ^ {zjrY inclusive. This is a rigorous result, since the coefficients 7(a)(f))(c) 
under discussion have no derivatives of the basis components with respect to z, 
and therefore the possible unaccounted for corrections due to terms ^ (z/r)^ in 
Cz and cannot contribute. Second, direct calculation similarly shows that 
l(t)(z){z) = 0 through the order -- (z/r)^. 

The calculation of all non-zero coefficients of the second type is much more 
cumbersome task. But as we will see below, the only coefficient of this type 
that is needed has the form 


z H 

i(z)(t)(t) — 

We note that all connection coefficients of the type 7(a) (t)(t) vanish in the 
equatorial plane z = 0. This is consistent with the requirement that the four- 
velocity 17^“^ = (1, 0,0,0) must satisfy the geodesic equation at z = 0: 

-^ = C/'’e(,) ([/(“)) + = 7(a)(t)(t) = 0- (1.58) 

1.2.3 Relativistic hydrodynamic equations 

Everywhere below, we only use the tetrad components of vectors, tensors and 
covariant derivatives. Therefore, starting from this Section, we will substitute 
the tetrad notation by the standard one, which is familiar when using the co¬ 
ordinate basis. This means that from now on we do not put tetrad indices 
in parentheses and denote them by Latin letters *,j, /cEl; also, we denote the 
connection coefficients by T. 

The stress-energy tensor of a viscous fluid with energy flux has the form (see, 
e.g., paragraph 4.3 in |13j or paragraph 22.3 in M) 

= (p -F e -f p)WU'^ - -h 2?7cr*'= -H - C/'=g^ (1.59) 

^If one of the symbols (p, r, 2 , appears among the indices, it means that the corresponding 
index takes this value. 
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where p, e, p, p, C, are the rest-energy density, internal energy density, pressure 
and two viscosity coefficients, respectively, as measured in the local comoving 
fluid volume; q is the energy flux inside the fluid as measured by the local 
comoving observer. 

The shear tensor is 

^ik ^ 1 ( 1 - 60 ) 

with the projection operator 

pik ^ ^^k _ JJ^JJk^ 

The divergence of four-velocity is 

& = (1.62) 

The relativistic Euler equation is written as 

P^sT«^fc=0. (1.63) 

The energy conservation law has the form 

C/«T*'=;fc=0, (1.64) 

the rest-energy conservation law reads 

(pt/");fe = 0. (1.65) 

The covariant derivative in a non-coordinate basis is 

while the divergence of a rank-2 contravariant tensor is 

The energy flux vector and the the shear tensor (the deformation tensor free 
from pure scaling) are purely space-like objects: 

Uiq^ = 0 , = 0 , a\ = 0 . ( 1 . 66 ) 

1.3 Construction of the standard accretion disk 
model 

1.3.1 Basic assumptions and the vertical balance equation 

Thus, we consider a disk from the standpoint of local observers rotating around 
a black hole near its equatorial plane with a relativistic Keplerian velocity. 
Before writing the dynamic equations in the projection onto tetrad (11.49111.5^ . 
we discuss basic assumptions of the model and their consequences. In addition to 
obvious assumptions about axial symmetry and stationarity of the flow (meaning 
that the derivatives dt and are zero) the main assumption, which we have 
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already used, is the small disk thickness, <5 = h{r)jr <C 1, where h{r) is the 
characteristic height of the disk along the z-axis (more precisely, the disk half¬ 
thickness) . 

The disk symmetry with respect to the plane z = 0 implies that U*,U^, , q*, q^, g’', p,p, r], (, e 

are even functions of z, and and q^ are odd functions of z. 

We also assume that the characteristic scale of variations of these quantities 
in the radial direction is much larger than that in the vertical direction, that is, 
their ratio is greater than ^ 

Next, kinematic arguments suggest that 

- 6Ur (1.67) 

If the energy flux determined by vector q is proportional to the internal energy 
gradient e, then, for the local comoving observer, q*^^ = 0 and qf^^Tq^^^ 

Taking (ll.67|) into account, implies that the projection of q onto the four-velocity 
of circular equatorial motion is also small, i.e. of the order of ~ ^Qfoc- From 
the standard Lorentz transformations, we obtain that q*,q‘^,q^ ^ dq’^, i.e. the 
energy flux relative to the tetrad should be directed mostly normally to the disk 
plane. 

Now, taking all the above into account, we consider the projection of the 
relativistic analog of Euler equation (11.631) onto the ort in more detail: 


+ = 0. (1.68) 

Using the symmetry of physical quantities discussed above and symmetry 
properties of tetrad orts and connection coefficients (which become odd func¬ 
tions of z if they have at least one index z), discussed in Section fl.2.21 it is easy 
to check that equation (11.6811 is indeed an odd function of z. Further, we see 
that the first term in (11.6811 separately yields the term dzP and other terms con¬ 
taining p are smaller due to the smallness of U^. All other terms together can 
always be written as ^ zp/(r)(l -|- g{r,z)) with the function g{r,z) ^ 0{S^). 

Thus, we arrive at the important conclusion that necessarily 


1 dp 
p dz 


5 < 1 . 


(1.69) 


This means that in a thin disk the variables p, drP ^ i.e. these vari¬ 
ables are small relative to the dominant action of the gravitational force in this 
direction. Therefore, particles of the disk must move in trajectories close to 
geodesic ones. Clearly, in a steady-state and axially symmetric flow, this can be 
realized only in two cases: when the matter moves almost radially towards the 
gravitating center (and the specific angular momentum in the disk is everywhere 
close to zero) or when the matter moves in almost circular orbits (and the spe¬ 
cific angular momentum, oppositely, is maximal). We note that both cases are 
consistent with the general assumptions discussed above and the result (11.691) . 

^We note that the assumption about the velocity components in the disk plane, (7’", U‘^, 
should also be made that in vertical direction they can substantially change only on scales 
~ r, otherwise, terms in the shear tensor could arise that strongly dynamically contribute to 
the vertical balance condition, which would lead to a disk totally different from the basic case 
of interest here. 
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However, in the last case, the strict vertical hydrostatic equilibrium holds in the 
disk in the first order in (5; in other words, (11.681) can be rewritten in the form 

zf{r){l + 6 ^...). (1.70) 

p az 

When the flow is almost radial, the corrections in the parentheses in (11.7011 
are not small, and their value is determined by the contribution from the pre¬ 
vailing radial motion, when, due to the change in the disk thickness at each 
radius, the particles are accelerated in the z direction. 

Thus, the standard disk model includes one more independent assumption 
on the closeness of the fluid particle trajectories to equatorial circular orbits 
around the central black hole. Therefore, we will additionally suppose that in 
our reference frame ^ sU* with s <C 1 and later we can see how this 

second small parameter is related to 5. 

Consequently, we write equations first not only in the leading order in 5 but 
also by assuming s = 0 , i.e. that the flow moves along geodesic orbits and = 
(1,0,0, 0). Wherever needed, we then additionally evaluate the contribution 
from the terms in the leading order in s. 


Deformation of the velocity field 

We first find the non-zero components of the shear tensor in the leading order. 
First, the velocity divergence vanishes: 


0 = w,, = = r^oj = 0 


(1.71) 


Next, we have 

-ri,, 

and, in view of the symmetry in i and k, we see that the only non-zero compo¬ 
nents of the shear tensor are are 


r<j> 


= -l{r‘,.. + r„) = 1(1^+ .=/») 


3 D 

4 rP^C’ 


(1.72) 


a- = -ir^^ = 0(^). 


(1.73) 


Equation of hydrostatic equilibrium 

Substituting {7® = (1, 0, 0,0) in (11.681) and taking the smallness (due to the low 
sound velocity in the flow) of several non-zero terms containing 77 and compo¬ 
nents of q into account, we obtain 


dp 

dz 


pn 


z H 
k 3 ■ 


(1.74) 


Radial direction 

The radial projection of the relativistic Euler equation for s = 0 reads 

= 0, (1.75) 
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and excluding terms ^ containing the connection coefficients and components 
of q, we have only one non-zero term of the order of which has the forir(^ 

Clearly, this term should be balanced by the leading terms ~ s. Evidently, 
the contribution from 

[pU^U%k 

should be considered first, and here it can be only due to terms containing one of 
the connection coefficients of zeroth order in z and the time velocity component. 
There is only one such term: 2T1^U*U‘^ = 

Hence, we reach an important conclusion that s ^ 6^, i.e. the velocity 
components in the disk plane are 




(1.76) 


which is used when determining the force balance in the azimuthal direction. 


1.3.2 Azimuthal direction 


We consider the last projection of the relativistic Euler equation, its component 
along the azimuthal ort. Let us proceed in the same way as above and first 
write terms that are present in the case s = 0. Again, we take t/* = (1,0, 0,0) 
and see that 

[{p + e+p)U'^U%k=0, 


because = 0 through the order ^ (see the discussion at the end of Section 
ll.2.2p . Next, the term with pressure is absent by virtue of the axial symmetry, 
and terms with q® cannot contribute to the order higher than ~ <5"^. 

It thus remains to consider the contribution 


[2r,a*% = + {2r,a^^),, + + pO{S^) = 

= +>)0(i"). (1.77) 

Here, we are also dealing with terms of the second order in S^; therefore, it 
is necessary to hnd the leading contribution from terms ~ s. Again, we consider 
only the prevailing part due to ideal fluid term: 

- U^U,)[pWU%k. 

The second part, which is proportional to can be neglected because the 
term in square brackets cannot contribute to the zeroth order in 5, since there 
are no connection coefficients of the form TJ^ ~ (5°, as was discussed at the end 
of Section fl.2.21 

®The order of components can be estimated as follows. In the stationary case, the 
divergence of the energy flux must be of the order of the power generated due to viscous 
dissipation, which is, in turn, proportional to some scalar characterizing the degree of the 
velocity shear and the viscosity coefficient rj. In our case, the viscosity coefficient 77 < phcs 
(5^. The divergence is mainly due to the term dzQ^. This immediately implies that ~ <5^ 
and 5^. 
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As a result, we obtain 


[pu*u\k = pvtu^u^ + prffct/^c7' = p{Tt + rfjc/'- = 


U'^ flH 

1^2 C 


U’' 

^2^C’ 


(1.78) 


where 



8 a 

j,3/2 


We now introduce the notation 


3a^ 

J.2 ■ 






(1.79) 


(1.80) 


where is the vertically integrated density of the flux of the (/)-coniponent of 
momentum in the radial direction. Then, by integrating ()1.77ll and (11.7811 over 
the disk thickness and combining them in one equation, we have 


^ ^ A _ E ^ 

dr rD \ r J 2r3/2 CD^/'^ ’ 


(1.81) 


where the contribution from vanishes due to its being an odd function of 
z, and we neglect the dependence of on z, which gives rise to a higher-order 
correction (see footnote 2). In formula (ll.8ip we have introduced the surface 
density of the disk 



(1.82) 


Important equation (11.811) with known boundary conditions at the inner disk 
radius allows us to calculate the profile T,y{r) for disk provided that the radial 
velocity distribution is known. 

We note that equation for T^, can also be derived from the angular momen¬ 
tum conservation law, which was used in the original paper [U (see equations 
(5.6.3)-(5.6.6) therein). 


1.3.3 Rest energy conservation law. Radial momentum 
transfer 

To solve equation (I1.81L the radial velocity profile should be specified. It can 
be obtained from the rest energy conservation law (11.651) : 

eripUn+e.{pU^) + ri,pU^ = 0. (1.83) 

Clearly, the substitution {7* = (1, 0,0,0) does not yield non-zero terms up to 
the order ^ (see the discussion at the end of Section 11.2.21) . In our reference 
frame, this fact can be easily understood: the circular axially symmetric motion 
corresponds to zero velocity divergence. It is straightforward to check that the 
following terms ^ s will appear in the continuity equation: 

D^^^pU^),r + (pun.z - pU^ = 0, (1.84) 

r 
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where the last term arises due to the contribution from and similar 

terms with other velocity components, even if they appear, have an order higher 
than ^ 5'^. 

After integrating over z, the contribution from the second term in (11.8411 
vanishes because p —>■ 0 far from the equatorial disk plane, and we obtain 

= 0. (1.85) 


The combination whose derivative is found in H1.85|) is a constant, which 
must be identified with the radial flux of matter. After additional integration 
over (j) we obtain that 

=-M, (1.86) 

where M > 0 is the rate of the matter inflow into the disk at infinity, i.e. the 
mass accretion rate. 

After substituting (11.8611 in (11.8111 . we finally obtain 


dT, 

dr 


PiTu + P 2 — 0, 


(1.87) 


where 



_ M E 
^ drr r^/'^CD 

The solution to ()1.87ll with the boundary condition T\r^^ = 0 is written in 
the form 


T. = 


P 2 {x)F{x) dx, 


F{r) 

F{r) = exp Pi{x)d3 


The integral (11.8911 is elementary, and as a result we obtain 

M r E 


= 


Itt r'^D 


r^/^C 


dr. 


( 1 . 88 ) 

(1.89) 


(1.90) 


1.3.4 Energy balance 

Here, we consider equation (11.6411 . As above, let us set 17* = (1, 0, 0,0) and find 
terms of the leading order in (5. As in the case of the azimuthal projection of 
the relativistic Euler equation, ‘ideal’ terms [(p + e + p)U*U^].k and prp^.^k 4o 
not contribute here. From the shear term, we have 


[2p(T^%k = = 2p[{Tt^r + + 0 (^')] = + 0 (^')]. 

Terms with g* contribute due to rapid change in the energy flux component 
normal to the disk with z: 

(C/V);fe = ^ + 0((54)- 
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Summing all terms, we obtain from the energy balance equation 


dq^ 

dz 




(L91) 


whence, after integrating over the disk thickness, we derive an important relation 


Q 


3 D 




(1.92) 


where Q = q^{z = h) is the vertical energy flux escaping from the disk. After 
specifying Q, we can calculate the radial profile of the effective temperature of 
the disk surface, because by definition Q = This is the universal result 

of the standard accretion disk theory: Teff does not depend on the specific 
nature of the dissipation of the kinetic energy of matter or on the mechanism 
of thermal energy transfer toward the disk surface, and is proportional to the 
value of M, times some universal known function of r. 

Thus, we have obtained the explicit form of the viscous stress integrated 
over the disk thickness, Tb, and the explicit form of the radiation energy flux 
from its surface, Q. At the same time, we know only the combination EC/’', 
and not each of these variables separately. In addition, we should determine 
the disk half-thickness profile, /i(r), and the temperature, pressure and density 
distributions, T{r,z), p{r,z) and p(r,z), inside it. To do this, the vertical disk 
structure should be calculated. 


1.3.5 Energy transfer equation and the vertical disk struc¬ 
ture 

The vertical disk structure is defined by three equations. Two of them have 
already been obtained above: the vertical hydrostatic balance equation (jl.74l) 
and thermal energy generation equation (11.911) . 

The remaining equation is the transfer equation for energy dissipating in the 
disk. In the simplest case, the energy transfer is due to the photon diffusion 
in heated matter. Strictly speaking, we should writ a relativistic analog of the 
radiation heat conductivity equation, which is a variant of the kinetic Boltz¬ 
mann equation for photons when their mean free path length is much smaller 
than the characteristic spatial length of the problem. Boltzmann equation is 
relativistically generalized in section 2.6 in [I]. The standard transition to the 
diffusion approximation yields the following equation (see expression 2.6.43 in 

my- 

q^ = ~bT^P^\ek{T) + UfcT), (1.93) 

Kp 3 

where R is the Rosseland mean opacity of matter, T is the temperature, b is the 
radiation constant and ak = Uk-jU^ is the four-acceleration. The discussion of 
equation (ll.93|) can be also found on p. 165 of [13]. 

As regards (11.9311 . we first note that the four-acceleration never exceeds the 
order ^ 5^, because the four-velocity itself differs from the geodesic value (free 
circular equatorial motion) only in the second order in S. In contrast, the 
derivative in the first term in parentheses in the right-hand side of (11.931) for 
k = z raises the order in S, since T, as well as e, , vary significantly across the 
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disk thickness. As a result, as already discussed in Section fl.3.11 we see that q 
is the leading component of vector q and is determined by the equation 


1 d{bT^) 

Skp dz 


(1.94) 


which is identical to the Newtonian form for a thin disk. 

Equations (11.7411 . (I1.91|) and (11.9411 must be supplemented with the equation 
of state of matter 

PiP,T), 


the opacity law 


k{p,T), 


and the explicit form of 

7?(p,r),or T;^p,r) 

depending on the type of parametrization of the turbulent viscosity in the disk. 

In addition, it is necessary to set boundary conditions at the integration in¬ 
terval z G [0, ft.]. In the simplest case, we assume that the disk has no atmosphere 
and 

p\z=h = Tl^^h = 0 . 


Furthermore, the energy flux vanishes in the disk equatorial plane: 


9*|z=o — 0. 


Finally, we denote 

ph 

2 / Tl'^dz = T^. 

Jo 

Note that the above equations and boundary conditions for the vertical disk 
structure automatically guarantee the validity of equations (I1.86L (|1.9()ll and 
(ll.92|l for the radial disk structure. 

After calculating the vertical structure, we can specify the surface density 
distribution using (11.8211 and then t/’’ using (11.8611 . 


1.3.6 Parametrization of turbulent viscosity and the ex¬ 
plicit disk structure 

Estimates carried out in m and [T] according to the algorithm described in 
Section fl. 3. 51 show that at sufficiently high accretion rate M, which is the free 
parameter of the problem, the radiation energy becomes dominant in the inner 
parts of the disk. The estimate of the threshold value of M can be found, for 
example, in [T^ (see formula 2.18 therein). It turns out that the disk thickness 
far away from its inner radius is independent on r and for M of the order of and 
above the critical value, Mcr (when the disk luminosity reaches the Eddington 
value in the inner parts of the disk), <5 > 1, corresponding to the spherization 
of the flow (see expression 7.1 and its discussion in m)- In addition, later 
studies showed that the radiation-dominated region is thermally unstable [16] 
and convectively unstable m 

This means that for the correct description of the inner parts of accretion 
disks at high accretion rates, when S increases, terms of higher order in S should 


23 























be taken into account. These include the radial pressure gradient ^ <5^ in the 
radial force balance and the advection term, U^TdS/dr ^ (5^, which arises in 
the energy balance and accounts for the radial heat transfer. The latter, in 
fact, implies that the heat diffusion time in the vertical direction is comparable 
to its radial advection due to radial transfer of matter. In other words, the 
main property of the standard accretion disk model considered here is violated: 
the local energy balance in the disk, when the heat generated due to turbulent 
energy dissipation is locally released from the disk surface. It was found that 
the account for the new terms also allows one to correctly describe the region 
near Tms, where in the standard model C/’’ —oo, and to construct a stationary 
solution with i5 < 1 for M of the order and above Mcr that is stable under 
thermal perturbations (so-called ‘slim-disks’, see [TB] and m and their citations 
list, and, e.g., [lO])- Later, these results were confirmed by numerical simulations 
(see, e.g., m and m)- We add that the transition from a standard disk to a 
slim disk with increasing M in the relativistic model around a rotating black 
hole should occur even earlier due to a higher accretion efficiency (which is, in 
turn, due to both decreasing rms and the additional angular momentum loss 
from the disk surface by radiation). 

Now, assuming that M <C Mcr, let us calculate the disk vertical profile, 
which is to be useful in the next part of the paper, in the simplest case where 
the pressure is mainly determined by fully ionized hydrogen plasma, i.e. 

p = 2pkT/mp, (1.95) 

where nip is the mass of a proton, is the Boltzmann constant, and the opacity 
is determined by Thomson scattering, k = kt = QAcw? jg. 

Let us also assume that the kinematic viscosity v is independent of z and 
can be parametrized in the form 


v = acsh, (1.96) 

where 0 < a < 1 is the Shakura parameter determining the turbulent viscosity 
in the disk (see [13] and m), and Cs is the speed of sound in the equatorial 
disk plane. Here, due to (I1.95L 


cl = 2kTc/nip, 


(1.97) 


where Tc = T(z = 0). 
Equation (11.9411 yields 


dzq^p=-—bT^ 

6Kt 


Zkt 


-bTt 


On the other hand. 


rn ^ 

dzq^p = CqF / pdz = -Cq'EF, 


where Cq is some correction factor of the order of unity corresponding to the 
difference between the escaping radiation flux, Q, and its mean value along the 
disk thickness. As a result, we get 
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To = ^ 

Next, for simplicity we assume that the entropy is constant along z and, 
dividing the left-hand side of (jl.74l) by p, we introduce the enthalpy, dw = dp/ 
integrate (11.7411 over z and obtain the central value of w, Wc = w(z = 0): 

_ Jf_H 
r3 C “ 2r3 C ’ 

Hence, using that Wc = nc/., where n is the polytrope index, we obtain 



c? = 


/l2 H 

2nr3 C 


(1.99) 


Finally, due to definition (11.801) . parametrization (11.961) and equation (11.901) 
we find 


^ 3 74 MY 


2r3/2C' 27rr3/2Z4’ 

where in the second equality we introduce the new variable 




dr. 


( 1 . 100 ) 


( 1 . 101 ) 


which in the Newtonian limit, far away from the inner edge of the disk, tends 
to unity. 

Equations (I1.92L (11.9711 . (11.9811 and (11.9911 are sufficient to exclude all un¬ 
knowns except S and free parameters M and a from (11.10011 . We thus obtain 
the following surface density profile E: 

E = (1.102) 


where the dimensional constant Eg combines all relevant physical constants and 
numerical coefficients. Its explicit form and numerical value (which depends on 
the black hole mass to which we normalize all quantities) can be found by the 
reader. 

Now, using formulas (I1.99L (11.9711 . (11.981) and (11.10211 . it is possible to derive 
the profile h{r). The resulting disk aspect ratio reads 5{r) = h{r)/r: 

S{r) = sy/20(j9/20j^-i/5jj-9/20Yi/5^ (1.103) 

where <5* is a constant that determines the characteristic disk thickness S. 
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Chapter 2 

Relativistic twisted 
accretion disk 


2.1 Introductory remarks 


In Chapter 1 we described a flat disk in the equatorial disk around a rotating 
black hole. Its axially symmetric structure was evident and consistent with the 
symmetry of space near the black hole. If now we relax the main assumption 
that the flow of matter at all distances coincides with the equatorial plane, the 
question arises: what can the dynamics of this more complicated flow, both 
stationary and non-stationary, be? Is this flow similar to a disk in any way? 
For thin disks considered here, the answer to this question proves to be positive 
under some restrictions. 

The main reason for the deformation of (for example, initially flat) disk 
is that the black hole spin gives rise to an additional off-center gravitational 
interaction with the gas elements of the flow. It can be shown that far away 
from the event horizon but close to the equatorial plane of the black hole this 
interaction is presented by an axially symmetric field of force directed to the 
black hole spin axis in the planes parallel to the equatorial one (see [53] chapter 
3, paragraph A). 

This force is called gravitomagnetic force and is given in this case by the 
expression 


Fqm 


dofl d 

r2 


( 2 . 1 ) 


where fl is the Keplerian frequency and d/dr is the radial coordinate ort of 
the cylinder reference frame. Clearly, this external force can change the proper 
angular momentum of the disk elements (and hence deform the disk) moving 
outside the equatorial plane of the black hole. Here only the projection of the 
gravitomagnetic force onto the angular momentum direction matters, which is 
proportional to the sine of the angle between the angular momentum vector and 
the black hole spin axis. As we will see shortly, the restriction that allows us 
to treat the new configuration as a disk (both stationary and non-stationary) 
requires that the gravitomagnetic force be smaller than the central gravitational 
attraction force, i.e. requires the smallness of the parameter a <C 1. In addition, 
one more restriction can be formulated that the non-complanarity of the disk 
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with the equatorial plane of the black hole, as well as the degree of its deviation 
from the planar form (i.e. twist, warp) should not exceed some small values for 
the disk to be hydrodynamically stable (see [33] paragraph 7 and [S] paragraph 
4.2.4). 

Let us split a thin planar disk into rings of small widths. In each ring, the 
motion of gas elements is mainly due to gravitational attraction force from the 
central body. The characteristic time of this motion is ~ 12“^. In addition, 
td determines the time it takes for the disk to restore the hydrodynamic equi¬ 
librium across the ring, since the disk aspect ratio (the ratio of the the disk 
thickness to the radial distance) is of the order of the ratio of the sound velocity 
to the orbital velocity. This conclusion can be also arrived at by noticing that 
vertical pressure gradient is times smaller than the unit mass gas element 
acceleration, i.e. exactly as small as the ratio of the radial size of the ring to its 
vertical scale. Thus, we can conclude that if other forces acting on a given ring 
from the adjacent rings or from the black hole lead to the dynamics with the 
characteristic time tey much greater than the dynamical one, t^v ^ td^ the hy¬ 
drostatic equilibrium is maintained in the ring, in other words, the ring remains 
flat, and the entire flow preserves the disk-like form. This is undoubtedly so in 
a flat disk, because in this case equally oriented rings interact by the viscous 
force acting in the azimuthal direction and the angular momentum changes due 
to inflow and outflow of the matter accreting through the ring, with the both 
processes occurring on the diffusion time scale, ty ^ td- 

Now let the disk be tilted with respect to the equatorial plane of the black 
hole by a small angle /3 <C I. If in a flat disk the gravitomagnetic force con¬ 
tributes only to the modulus of acceleration of gas elements moving in circular 
orbits, but now, due to a non-zero projection of this force (oc /3) onto the an¬ 
gular momentum of gas elements, this force makes the orbits to process around 
the black hole spin axis. For free particles, this effect is described in detail 
in the second part of the next Section in terms of the difference between the 
frequencies of circular and vertical motions. We also show that the precession 
frequency is much smaller than the circular frequency for a <C 1 (see formula 
(I2.12|l l. which is equivalent to the condition tev ^ td for whole rings composed 
of gas elements. 

Equation (12.121) suggests that the precession of the rings is differential, i.e. 
depends on the distance to the center. As a result, the relative orientation of 
initially coaxial rings changes and the disk is no longer flat. However, we keep 
in mind that under the condition tev ^ td each of the rings behaves ’rigidly’ in 
its vertical direction, which is now also a function of r. The new configuration 
is similar to a twisted (or warped) disk, i.e. a flow symmetric relative to some 
(now not planar) surface, which can be called the equatorial surface of the 
twisted disk. Here, the cross-section of the equatorial surface by a plane passing 
through the center is a circle - the instantaneous shape of orbits of gas elements 
rotating with a given radial distance r. The disk turns into a set of rings tilted 
to the black hole equatorial plane by a constant angle /3 but with depending 
on r node lines (the line formed by the intersection of the ring planes with the 
black hole equatorial plane). The node line is now determined by the position 
angle 7 (r) measured in the equatorial plane in the positive direction from a fixed 
direction to the ascending node of a given ring. The key point here is that the 
pressure gradient in the twisted disk, directed (as in any thin disk in general) 
almost normal to its warped surface, is not normal to the planes of the rings 
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composing the disk. Therefore, we conclude that the pressure gradient acquires 
two projections. The main projection is coaxial with the rotational axis of each 
ring. Let us conventionally denote it as (Vp)j, where ^ is the distance from the 
equatorial surface of the twisted disk measured along the direction of rotation of 
the ring reduces to z in the case of a flat disk). We note from the beginning 
that (Vp)^ oc ^ due to the hydrostatic equilibrium across the ring. The second 
projection of the pressure gradient, conventionally denoted as (Vp)r, lies in the 
ring’s plane along the radial direction connecting the disk center and a given 
gas element of the ring. The ratio of these two projections is a small value 
proportional to the rate of change of orientations of rings in the disk, which, 
in turn, depends on the radial direction chosen in the given ring’s plane. From 
purely geometrical considerations, we rigorously show in what follows that for a 
disk with /3 = const the ratio (Vp)r./(Vp )5 oc Pdj/dr cos tp, where tp is the angle 
measured in the azimuthal direction for the given ring from its ascending node 
to the given gas element. Note that the normal to the twisted disk surface is 
orthogonal to the ring’s plane only in two diametrically opposite points - where 
the ring’s plane intersects the planes of the adjacent rings. At /3 = const, these 
points are characterized hy ip = ±tt/2. At the same time, in other pair of points 
with Ip = 0,7: the value of {Vp)r reaches positive and negative maxima. 

Thus, in a flat disk, the dynamics in the radial direction is controlled in the 
leading order in J <C 1 by the gravitation force and the corrections ^ S'^ are 
neglected, whereas in a twisted disk the radial projection of the pressure gradient 
starts additionally contributing to the radial balance. This addition, on the 
one hand, depends on the degree of the twist, and on the other hand, increases 
proportionally to the distance from the equatorial disk surface, Next, because 
it also depends harmonically on the azimuthal direction, the gas elements (for 
^ ^ 0) are subjected to periodic disturbance by this force with the orbital period, 
and their orbits become ellipses with small eccentricity. As is well known, the 
eigenfrequency of the small oscillations of free particles in eccentric orbits is 
equal to the epicyclic frequency, k. Because the pressure gradient projection 
considered here excites exactly such oscillations, the radial profile of the epicyclic 
frequency, K(r), is an important characteristic that determines the shape of both 
stationary and non-stationary twisted configurations. In the next Section, we 
derive the required relativistic profile K(r) for equatorial circular orbits in the 
Kerr metric (see equation (12.81) 1. Note from the beginning that in the special 
case of Newtonian gravitation k = 17, and hence the action of the external 
exciting force on gas elements with the same frequency results in a resonance: 
the amplitude of the perturbed motion, characterized by perturbation of the 
orbital velocity, v, must increase without a bound. This growth, however, is 
always limited by turbulent viscosity in the disk. Indeed, since the exciting force 
amplitude oc so is the amplitude of v. But this would mean the presence of 
the vertical velocity shear, d^v, in each ring. Together with the vertical density 
gradient (and hence the vertical gradient of the dynamic viscosity) in the disk, 
this gives rise to a volume viscous force that damps the driving of individual 
layers of each disk’s ring by the resonance force. Note that near the black 
hole, where the frequency n deviates from 17, the amplitude v remains bounded 
even in the absence of viscous forces. This allows the existence of stationary 
twisted disks with low viscosity around the black holes, in which /3{r) takes an 
oscillatory form (see [25] 1. 

Thus, we see that the twist of the disk caused by the gravitomagnetic force 
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necessarily results in a perturbation of the circular motion of gas elements in 
disk’s rings. The velocity field of this perturbation, v, depends on r (in addition 
to its being proportional to oc as explained above) and is determined by the 
current shape of the disk. By virtue of the continuity of the flow, this gives rise 
to density inhomogeneities outside the disk equatorial surface, pi oc Because 
(Vp)r oc cos these inhomogeneities take the opposite signs in the diametrically 
opposite points of any given ring. But this implies that the ring is subjected to 
the total torque of the central gravitational force acting on the density excesses 
out of the equatorial plane of the ring (i.e. outside ^ = 0). We denote this 
torque as Tg. Because the disk is thin and the gravitational acceleration along 
ring’s axis is itself oc the corresponding component of the gravitational force, 
and Tg as well, are quadratic in In addition, we remind that the torque Tg 
is proportional to the small warp magnitude, Tg oc pd'y/dr. Thus, we arrive 
at the conclusion that the dynamics of the twisted disk rings is controlled by 
Tg, together with the torque due to the gravitomagnetic force discussed earlier 
in this introductory Section. Note that in the case /3 = const considered here, 
(Vp)r and, correspondingly, pi take the maximum absolute values (but with the 
opposite signs) at tjj = 0, tt, i.e. at the node line of each ring[3. But this implies 
that Tg lies in the plane made by the angular momentum of each ring and 
the black hole spin axis. By virtue of the symmetry of the problem, the total 
contribution to Tg from other azimuths does not alter its direction. Therefore, 
immediately after the gravitomagnetic force turns an imaginary tilted planar 
disk into a twisted configuration with j3 = const, the gravitational force acting 
on the asymmetrically located matter of the disk relative to the surface ^ = 0 
tends to change the disk rings tilt angles: either to align the rings with the 
equatorial plane of the black hole or, conversely, to remove them from it. On 
the other hand, once /3 becomes dependent on r, the maxima of absolute values 
of {Vp)r are shifted from the node line of each ring to some new tp, which gives 
rise to a component in Tg that also contributes to the precession motion of the 
disk rings, as the gravitomagnetic torque does. 

The dynamics of twisted disks sketched above is complicated by the presence 
of non-zero viscosity in the disk. First of all, each ring of the disk is subjected 
to the action of the viscous force arising due to the difference between the 
direction of the tangential velocity of the ring and that of the adjacent rings. 
This difference is maximal in the directions where the ring planes intersect, 
i.e. exactly where (Vp)i. vanishes. In the above example of configuration with 
/? = const this corresponds to ip = ±7r/2, i.e. perpendicular to the node line of 
the rings. The viscous force, being proportional to the difference in tangential 
velocities, is directed at these points perpendicular to the ring plane and has 
different signs on different sides from the node line. Therefore, the corresponding 
torque, Ti, is perpendicular to the plane made by the ring’s angular momentum 
and the black hole spin axis. In other words, the viscous interaction between 
the disk rings leads only to their precession around the black hole spin. Note 
also that the viscous torque Ty oc jdd'y/dr, which appears due to the difference 
between the tangential velocities of adjacent rings, and Ty oc due to the 
viscosity coefficient. It is important to note that as soon as the profile /3(r) is 
formed due to the gravitational torque Tg, the viscous torque Ty starts causing 

^To make the description as rigorous as possible, it is important also to add that the noted 
coincidence of azimuthal location of maxima of (Vp)r- and pi occurs only when the action of 
viscosity on the gas elements of the ring is neglected. 
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alignment/misalignment of the ring with the equatorial plane of the black hole. 
This happens for the same reasons by which Tg also starts contributing to the 
precession motion, as discussed above. 

In addition to causing the appearance of T^, the viscosity in a twisted disk, 
like in a flat accretion disk, leads to the radial diffusion transfer of the angular 
momentum component parallel to the equatorial plane of the black hole (which 
is non zero exactly for a tilted/twisted disk) toward the disk center due to 
simple transport of the accreting matter, and toward its periphery due to the 
corresponding angular momentum outflow. In the case of a relativistic disk, an 
additional loss of this angular momentum component occurs due to the thermal 
energy outflow by radiation from the disk surface (see equation (C6) in [3). 

All forces participating in the dynamics of twisted disks appear in the so- 
called ‘twist’ equation - the principal equation of the twisted disk theory. This 
equation is derived and analyzed in the subsequent Sections. 

2.1.1 Weakly perturbed circular equatorial motion: epicyclic 
frequency and frequency of vertical oscillations 

In a twisted disk, the motion of matter outside the equatorial plane of the Kerr 
metric is assumed; this motion is not necessarily circular in the projection onto 
that plane. Therefore, we first analyze the properties of free particles moving 
in orbits slightly different from circular ones. 

We first assume that particles move exactly in the equatorial plane but 
in slightly non-circular orbits. The problem can be solved using relativistic 
hydrodynamic equations with zero pressure and by assuming that there is a 
small addition to the purely circular velocity. Then, instead of equations (11.631) 
and (|1.64L it is better to use the original equations in the form 

= 0 , ( 2 . 2 ) 

where in the considered case of free motion, = pWU^ and p = const. Under 
the last assumption, the velocity field, as follows from the rest-energy conser¬ 
vation law (|1.65|) . is divergence-free, and (12.211 is equivalent to the following 
equation: 

= 0. (2.3) 

We now single out from the four-velocity field a small addition to the main circu¬ 
lar equatorial motion and denote it as u*. The unperturbed motion corresponds 
to rest in the projection onto tetrad (11.49111.521) used to construct the flat accre¬ 
tion disk model, i.e. is given by the four-velocity Uq = {1,0,0,0}. Substituting 
the sum Uq + Vi in (12.31) . we obtain linear equations for small perturbations of 
the four-velocity, Vi, which is assumed to be a function of t only: 

v^-,kUS + = 0. (2.4) 

Taking into account that C/g.^. = T)^, for i = 1,2 we obtain the system of 
equations 

u’';* -h = 0, (2.5) 

uft + Tfy = + r-3/2 Q (2.6) 
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It follows that small perturbations of the four-velocity components in the equa¬ 
torial plane of the rotating black hole oscillate in time. For example, v'^ satisfies 
the equation 


(Pv’’ 

dt'^ 


2C 

r^B'^ 



(2.7) 


which implies that the square of the frequency of these oscillations, which is the 
epicyclic frequency by definition, has the form 


^2 = r-^B-^{2C -H)= r-3 (l + ^ 


-2 



r 


8a 

^ 3/2 



( 2 . 8 ) 


A somewhat different derivation of k can be found in the Appendix in [26) . 
It is important to note that (EHl) contains a derivative with respect to the 
coordinate time, and therefore the epicyclic frequency is determined by the 
clock of an infinitely remote observer, similarly to circular frequency (11.1611 
introduced above. By comparing equation (I1.17L which defines the location 
of the innermost stable circular equatorial orbit in the Kerr metric, r^sj with 
(El, we infer that K^{rms) = 0. For r < Tms the epicyclic frequency becomes 
imaginary, and equation (EB has exponentially growing solutions. It must be 
so because in this region the free circular motion around a rotating black hole 
becomes unstable. In Section fl. 1.31 this result was obtained from the analysis of 
the form of the effective centrifugal potential in which a test particle moves in 
an equatorial circular orbit. Nevertheless, we see that r^s can be determined 
alternatively from the calculation of the profile ^^(r) in the Kerr metric. 

It is well known that for Newtonian motion so-called Keplerian degeneration 
occurs when k = fl for a non-circular motion, which causes non-relativistic orbits 
to be closed. However, this symmetry is broken for relativistic free motion, and 
the epicyclic frequency k differs from fl already near a non-rotating (a = 0) 
black hole, where its square is 

(i - 7 ) =«“(i-;)<«’■ P-9) 


The difference between the epicyclic and circular frequencies results in the well 
known effect of the precession of an elliptical orbit. Far away from the horizon 
of a Schwarzschild black hole, i.e. for r ^ 1, the frequency of the orbit rotation, 
called the Einstein precession frequency, is flp ~ 

We now suppose that we rotate together with the test particle at some radius. 
When considering the problem in the projection onto tetrad (11.49111.,^ . this 
particle remains at rest. We now impart to the particle a small velocity in the 
direction perpendicular to the equatorial plane. Equation (11.741) of hydrostatic 
equilibrium for a flat disk implies that in our reference frame the particle, being 
in free motion, is subjected to acceleration that is proportional to z and tends 
to return the particle to z = 0. As a result, the test particle will harmonically 
oscillates with a frequency whose square is 


fit 


H 


( 2 . 10 ) 


where the superscript ‘1’ reminds us that the frequency is measured in the ref¬ 
erence frame comoving with the particle in its main circular equatorial motion. 
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To re-define this frequency as measured by the clock of an infinite observer, 
as has been done for both circular and epicyclic frequencies, the frequency Oj, 
must be divided by the time dilation factor (the difference between the proper 
time of the particle and the time at infinity), i.e. by the t—component of the 
four-velocity (11.141) . Thus, the square of the frequency of vertical oscillations is 


a,2 ^ = r 


„-3 


1 + 


-1 


' 3/2 / 


/ 4a 3a2\ 

I ~ j.3/2 ^2 ) ’ 


( 2 . 11 ) 


which coincides, for example, with the expression presented in m (see also 
[28]). Eq. (12.111) implies that around a non-rotating black hole ily = fl. This 
means that the vertical and circular motions have the same period, and the 
total motion of the particle is again the circular motion in a closed orbit whose 
plane, however, is now slightly tilted toward the initial equatorial plane. The 
situation changes for a yf 0, because for fly ^ the orbit is not closed any more, 
and the orbital plane starts precessing around the spin axis of the black hole. 
The frequency of the orbital precession is equal to the difference between the 
circular and vertical frequencies. For a slowly rotating black hole with o <C 1 
the precession frequency of a slightly tilted orbit is 


^LT = - ^ly 


j.3/2 






( 2 . 12 ) 


This is simply the angular velocity of the frame dragging by the rotating black 
hole (see equation da) in the limit a <C 1. The frequency ^lt is also referred 
to as the Lense-Thirring frequency. 

In the most general case, where the test particle deviates from circular mo¬ 
tion simultaneously in the vertical and horizontal directions, the particle’s mo¬ 
tion in space can be described by a slightly elliptical orbit, with both plane and 
apse line turning with an angular velocity proportional to the difference between 
the circular and vertical frequency and the difference between the circular and 
epicyclic frequency, respectively. For a <C 1, the precession of the orbital plane 
occurs on a timescale much longer than the dynamical time, ^ td, where 
Ilt discussion in the previous Section). 


2.2 Choice of the reference frame 

2.2.1 The metric 

Taking the general conclusions in Section 2.1 into account, we consider slowly 
rotating black holes, a 1. In this case, the linear expansion of the Kerr metric 
in the parameter a is sufficient. Then formula (HU takes the form 

ds^ = (1 — 2/R)dt^ ~ (1 ~ 2/R)~^dR^ — R?{d9'^ sin^Odcj)^) + d(j>dt. 

(2.13) 

Metric (I2.13P is identical to that of a non-rotating black hole written in the 
Schwarzschild coordinates, except for one non-diagonal term responsible for the 
Lense-Thirring precession. 

Our main purpose in this Section is to introduce the relativistic reference 
frame that follows the disk twist. The symmetry of the problem implies that the 
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equations of motion should have the simplest form in such a frame. As for a flat 
disk, it is convenient to use some orthonormal non-coordinate basis. For this 
basis to follow the disk shape, its two spatial orts should be tangent to the disk 
symmetry plane. At each spatial point we take the orts of the ’flat’ basis, which 
are determined, say, by the equatorial plane of the black hole, and turn them by 
the angles /3 and 7 defining the disk shape. This is done in the simplest way by 
using a Cartesian coordinate system with the axis z parallel to the black hole 
spin. However, we should first understand which four-dimensional basis (whose 
dual tetrad must transform the metric (I2.1dll into the Minkowski metric) in 
the flat-space limit would produce the spatial part described by the Cartesian 
reference frame. 

This can be done by changing the radial variable in (12.131) , namely, by passing 
from R to the so-called ’isotropic’ radial coordinate, Rj: 

Substituting (12.141) in (12.131) yields 

ds^ = (df- (1 + V {dRj+R^jd0^+Rjsin^ed^^)+4 — 

V + miJ V + 

(2.15)' 

where the second term represents the elementary spherical volume. Now, it is 
easy to transform to the Cartesian coordinates by the change {x = Rj cos (j) sin 0,y = 
Ri sin (j) sin 6,z = Ri cos 6}. With account for Rj sin^ 9d(j) = xdy — ydx we have 


ds^ = Kjdt^ + 2 aKiK 2 ,{xdy — ydx)dt — Kj{dx^ + dy^ + dz^), 
where 


Ki 


^ 2Ri 


K. 



2 1 



(2.16) 


(2.17) 


are functions of Rj = {x^ + only. 

Metric (12.161) generates the following dual basis 

e* = Kidt + aK^[xdy — ydx), = K 2 dx, = K 2 dy, = K 2 dz. (2.18) 


Note that basis (12.181) corresponds to the observers at rest in the Schwarzschild 
coordinates, because their world lines defined by the condition [/* = e'^/ds = 
{1,0,0, 0}, correspond to the equalities dx = dy = dz = 0. Their identical 
clocks are synchronized in such a way that in equal time intervals determined 
by the ort e‘, light travels an equal distance in any direction defined by the 
combination of the If the observers used the coordinate time, t, they 

would discover, for example, that the light signal in the azimuthal direction 
prograde with the black hole spin travels a larger distance than in the opposite 
(retrograde) direction. This follows from the frame-dragging effect of a rotating 
black hole and is equivalent to the well-known tilt of light cones in the azimuthal 


■^dt d(p, 
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direction. Finally, we note that another choice of the orthonormal basis is possi¬ 
ble in principle, which also compensates the space-dragging effect. Such a basis 
is called the frame of locally non-rotating observers, which is moving with the 
azimuthal angular velocity equal to (HU; mathematically, this corresponds to 
the correction of the azimuthal ort instead of the time one (see [29]). 

Below, we need to rotate the spatial part of (I2.18p . so as to obtain the dual 
twisted basis and then the original basis, which, as we recall, is needed to write 
down the projection of hydrodynamic equations. For this, let us first introduce 
the twisted cylindrical coordinates. 


2.2.2 Twisted coordinates 

We define the twisted cylindrical coordinates {t, r, tjj, such that the condition 
^ = 0 determines a coordinate surface coincident with the equatorial surface of a 
twisted disk. Here, r, r, tp and ^ are the new time variable and twisted analogs 
of the radial, azimuthal and vertical cylindrical coordinates, respectiveljH These 
coordinates were first introduced in [30] and [31]. At each fixed r = const, angle 
ip is measured in the positive direction from the ascending node of the circle 
^ = 0 crossing the equatorial plane of the black hole. The relation between 
{r, r, ip^} and {t, x, y z} can be obtained by a sequence of rotations at each 
radial distance by the angles p){r,T) and j(r,T). 

Let us take the radius vector with coordinates 


T 

rcostp 
r sin Ip ’ 


(2.19) 


where three spatial Cartesian coordinates are defined in a frame with the z-axis 
tilted by the angle /3(r, r) toward the black hole spin and the x-axis lying in 
the black hole equatorial plane and turned by the angle j(r) relative to some 
direction common for all r. 

Next, we consecutively rotate this frame by the angle pi{r,T) about its x- 
axis in the negative direction and then by the angle jir, t) about its z-axis in 
the negative direction. After these two rotations, this frame transforms into a 
’flat’ Cartesian frame common at all r with the xj/-plane coinciding with the 
equatorial plane of the black hole. Herewith, the new coordinates of the radius- 
vector are obtained by multiplying (I2.19P first by the matrix 


MP) 

and then by the matrix 

^2(7) 


10 0 0 

0 10 0 

0 0 cos /3 — sin /3 ’ 

0 0 sin /3 cos pi 


10 0 0 

0 cos 7 —sin 7 0 

0 sin 7 cos 7 0 

0 0 0 1 


^Here and hereafter, r denotes the twisted radial coordinate. 


( 2 . 20 ) 


( 2 . 21 ) 


34 













As a result, we obtain the following relation between the twisted cylindrical 
and ’flat’ Cartesian coordinates taken in the linear approximation in small /3: 

t = T 

X = r cos 7 cos '0 — sin 7 (r sin — ^/3) I'o oo'i 

y = r sin 7 cos '0 + COS 7 (rsin')/’— ^/3) ^ ’ 

z = r/3sin0 + ^. 

2.2.3 Tetrad transported by the twist-following observers 

We now pass from the ’flat’ basis (12.1811 to the twisted one by rotating its spatial 
orts by the twisting angles at each spatial point. First, we need to perform 
the rotation strictly opposite to what we did in the previous paragraph. This 
means that we should take basis (12.181) as a column and first multiply it by the 
matrix A 2 (— 7 ) and then by the matrix Ai(—/3). After that, because we wish to 
obtain the basis corresponding to the (twisted) cylindrical frame, it is necessary 
to additionally ’advance’ the three spatial orts by azimuthal angle 0 , what is 
achieved by additional multiplication of the basis by the matrix A 2 {—ip). 

As a result, we obtain the twisted dual basis that contains some linear com¬ 
binations of the ’flat’ coordinate orts, {dt, da:, dy, dz}. It remains to express 
it as linear combinations of coordinate orts of the twisted coordinate frame, 
{dr, dr, dip, d^}. For this, it suffices to take differentials of the coordinate trans¬ 
formation (given by ( 12 . 221 ) in the linear approximation in /3) and to substitute 
them in the twisted dual basis obtained after the rotations. It can be checked 
that in the approximation linear in /3 and a, we have 


e'’’ = {Ki — ar^K3d,pU)dT + a^K3d,p{Z—rW)dr+arK3{r — ^Z)d(f—arK3dipZd^, 

(2.23) 

=-^K 2 UdT + K 2 il-^W)dr, (2.24) 

= -^K2d^UdT - ^K2d^Wdr + rK2dip, (2.25) 

= rK2UdT + rK2Wdr + K2di, (2.26) 

where we introduce the new azimuthal variable Lp = ip + ^{r) and pass to partial 
derivatives with respect to the corresponding new coordinates. 

We also introduce new variables characterizing the disk geometry: 

'ki=/3cos7, 4'2=/3sin7 (2.27) 

and from now on use them instead of the angles /3 and 7 . Additionally, 

Z =/3sin0 ='hi sin(/j — 4>2 cos(/ 3 , U = Z, W = Z', (2.28) 


where partial derivatives with respect to r and r are denoted by the dot and 
the prime. 

It follows that for /3 = 7 = 0 and with the additional transition to the 
Cartesian coordinates, basis (12.23112.2^ is transformed into the ’flat’ basis (12.181) . 

As discussed above, observers transporting basis (12.181) are at rest in the 
Schwarzschild coordinates. On the contrary, observers corresponding to basis 
(12.23112.2S)l move in space by following the changing shape of the twisted disk 
(in the non-stationary dynamics). 
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As we have seen in Chapter 1, the original basis onto which hydrodynamic 
equations are projected is obtained by inverting the dual basis matrix. Using 
in the approximation linear in P and a, we have 


e, = ^ + pUdr + - rUd^'^ , (2.29) 

er = ^ (-ap^d^Zdr + (1 + pW)dr + - riua^) , (2.30) 

K 2 \ Ki r J 

^ ^ ~ ~ oC^rC/a^ ^ (r ~ ar^rUd^'^ , 

?2.31) 

^ (^ar^dy,Zdr + a^^ . (2.32) 

With the original and dual bases in hands, using the algorithm presented 
in Section 11.2.11 we can calculate the connection coefficients. This very cum¬ 
bersome but straightforward procedure yields the following non-zero connection 
coefficients in the linear approximation in /3 and a: 


K' 


Trrr — K 1 K 2 ’ 

TrrJ = —a^dipZ (1 — ^ -I- ^^4) , Tr<fir = —^Tr^p^ 


Trr^ = a^ {l-\{r-iZ)Ki), 


Trfr = 

r — ^^a jj — V 

-L ripr — UipKy x Tr(ft 

p _ (rK2)' _ „t K 3 a Tj 

i rcncn — ^ j, ^ 


p _ ^ 

4 tJt — J. K1K2 ’ 




nptp — 




r — ^-U-a lu 

4 npr — r K 2 ’ 


p _ U _p 

4 r{T — Ki 


p _ _ i ^2 

4 r^r - K 2 rK^ 


T- - K 


r<^{r — ^2 


K 1 K 2 

b(p{r = 

b<p{v> = -r^ ~ ar-g^do^U, 




(2.33) 

where K 4 = K^/Ki{Ki/K^)'. The other non-zero by^, as usual, can be ob¬ 
tained by taking their asymmetry in the first two indices into account. 

Thus, the basis (12.29112.3^ together with the connection coefficients (|2.3311 
are the sum of two parts: the main part that persists at ^ = 0 and a small 
additional part oc p. In what follows, we conventionally denote these parts 
“Bo” and “Bi”, respectively. 


36 














2.3 System of twist equations 

2.3.1 Projection of dynamical eqnations onto the twisted 
basis for a thin disk 

Separation of equations into two systems describing a flat disk and a 
twisted disk 

We take the relativistic hydrodynamic equations in the original form; 


T% = 0, (2.34) 

where the stress-energy tensor and its components are presented in Section ri.2.31 
Equations (12.341) should now be projected onto the twisted basis (I2.291I2.321) . 
To perform this, we assume that /3 <C 1. In other words, mathematically we 
consider the twist of the disk as a small perturbation to its ’ground’ state, i.e. to 
the model of a flat disk, also referred to as the background. It is important that 
the appearance of a twist gives rise to new terms in equations not only due to the 
bending of the basis, but also due to the appearance of additional perturbations 
of physical quantities themselves that enter the stress-energy tensor, including 
the density, pressure and four-velocity. 

For a twisted disk, instead of (12.341) we can write 


((ro*'= + + ((To*" + Ti*").,)i = 0 


ik 


ik 


ik \ 


(2.35) 


where Tq*" corresponds to the background state and Ti*" is a small Eulerian 
perturbation of the stress-energy tensor. The indices 0 and 1 that follow the 
notation of the covariant derivative mean that the derivative is taken in bases 
Bo and Bi, respectively. 

The action of the covariant derivative with index 0 on Tq*", evidently, yields 
0 , because these are equations for the background: 

(To*";fe)o = 0. (2.36) 


Then, in the approximation linear in /3, we find the twist equations: 

(Ti*".fe)o + (To*".fc)i =0. (2.37) 


We assume that in a twisted disk the four-velocity, pressure, rest-mass energy 
density, internal energy, viscosity coefficient and energy flux density, as defined 
in their standard sense (see Section 11.2.3p , are given by 

W = Uo+ v\ p = po + pi, p = pq + pi, e = Co + ei, 
r] = vo + m, q" = 90 + 91: 

respectively. Here, the indices 0 and 1 denote values related to the background 
and perturbations, respectively; u® are perturbations of the four-velocitjo 

Thus, Tg" is the stress-energy tensor that contains only unperturbed quan¬ 
tities in accordance with definition (I1.59P and its perturbation has the form 

T*" = wiU^U^ + woiv^U^ + C/^u") - Pi77*" + 2771(7^" + 2770 + 1 "- 

^To shorten the equations, we omit the term with the second viscosity (^: using the analysis 
be given below, it can be shown that this term does not contribute to the final equations in 
the leading order in the small parameters of the problem. 
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- U^<A - y'<lo - (2.38) 

where wq = Po + eo + Po is the background enthalpy and wi = pi + pi + ei is 
its perturbation. 

In addition, ctq^ is the shear tensor which contains only unperturbed quan¬ 
tities in accordance with definition (I1.60I1 . and is its perturbed part of the 
form 

+ (4)o^i - 

\[{Uo:,)oPt + (t^of,)oPn - liUoi,)oPl'^ + 

\moi,)iPt + iUofj)iPt] - l(.Uoi,)iPi\ (2.39) 

where Pq^ is the projection tensor that contains only unperturbed quantities 
in accordance with definition (ll.61|) . and its perturbation is written as Pl^ = 
-U^v'^ - U^v\ 

Everywhere below we omit the index 0 for the unperturbed quantities. In 
addition, the viscous part of the stress-energy tensor in the disk is marked with 
‘V” wherever necessary: = 2?7cr®^. 

Additional relations used to write the equations 

The relations given below are valid through terms of the order of oc 5^, which 
is sufficient for the theory of twisted disks in the leading order in the small 
parameter 5. In deriving these relations, such a simplification enables us to 
consider that in the background solution, only U'^ and t/''’ are non-zero, while 
t/’- oc 52 and t/’- can be temporarily set equal to zero. 

We first note that the following relation between the components U'^ and 
is used below: 

{U^f = {U^f + 1, (2.40) 

which follows from the expression for the norm of the four-velocity in an or¬ 
thonormal basis. Constraint (12.401) is also useful in the differential form: 

U'^dU^ = U^dU^. (2.41) 

Next, because the normalization of the four-velocity is also valid in the 
twisted disk, and the four-velocity perturbations are small, in the linear ap¬ 
proximation we have 

([/" + - {U^ + = {U^f + 2t/V - {U^f - 2U^v^ = 1, 

and hence, with account for (12.4011 . u* is “orthogonal” to t/b 

= U^v^. (2.42) 

Finally, from the condition that is space-like, we have 

and thus, in the basis Bg used in this Section, in the flat disk model, not only 
but also is non-zero in the order of <5 that is of interest for us here: 

T7 = ^T7- (2.43) 

Note that in basis (11.49111.5^ co-moving with the azimuthal motion, only the 
component was non-zero (see p.72p l. 
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Equation of free azimuthal motion 


The quantities corresponding to the background model and entering twist equa¬ 
tions (12.371) should be obtained separately from equations (I2.36|) . For this, it 
sufhces to use the results of Chapter 1 taking only the transition from basis 
(ll.49m.5^ to the basis Bq into account. 

Nevertheless, when deriving the twist equations, it is also necessary to use 
some of equations (12.361) written exactly in the basis Bq. We mean the r- and 
^-projections of these equations in the leading order in the small disk thickness 
which, as we know, describe its azimuthal rotation in the equatorial plane of 
the black hole and its vertical hydrostatic equilibrium. We emphasize that these 
relations are valid for both stationary and non-stationary accretion flow for any 
viscosity parametrization, as well as for any specific vertical and radial structure 
of the flow. Only the condition d <C 1 is important. 

At the first stage of deriving the twist equation we will need only the r- 
projection of (12.361) . Setting = pU^U^ we find that = 0 yields 


^{U^f + a—(2 - rK^WU^ - {jjvf = q. (2.44) 

Ki K 2 rK2 

Exactly this combination (12.441) is used in the derivation; however, it can be 
checked that together with (12.401) in the approximation linear in a it gives the 
solution 

= [rs - (1 - ar-'/2(rs - 3)-i) , (2.45) 

where we have transited to the Schwarzschild radial coordinate rg which is 
equivalent to r, that we used in Chapter 1 in the expression for (see formula 
(11.141) 1. It is easy to check that I/J = U'^ jrs, as must be the case with the 
transition from the coordinate basis to Bq taken into account. 


’Gauge’ condition of the twisted frame 

The principal kinematic constraint for the twisted reference frame requires a 
constant vertical position of fluid particles: 

f = 0, (2.46) 

which is provided by fast establishment of hydrostatic equilibrium across the 
disk compared to the dynamical time of the twist change, as discussed in Section 
2.1. However, as has been already noted in [32], an important point is that this 
does not mean that the projection of the four-velocity of the fluid onto is 
also zero, because our basis is non-coordinate and its orts are not tangent to 
the coordinate lines. 

By definition. 



Using (12.261) we have: 




dr 


dr 


= rK2U— +rK2W—] 


ds 


ds 
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we should substitute dr/ds and dr/ds in this relation in the zeroth order in /3, 
in other words, as values corresponding to the flat disk dynamics. Expressions 
for e^, and e’' at /3 = 0 give 


dr _ 1 / '^K dr _ dip _ U'^ 

ds Ki \ ^ ds) ' ds K2 ’ ds rK2 ’ 

where by definition W = e'^/ds. As a result, we obtain 

= rU^ K^U + rU^W, 

Ki 


(2.47) 


(2.48) 


where 


7,- I 1 

* = Ki IF 


In (I2.48|) the velocity components U'^ and C/’’ should be taken from the corre¬ 
sponding background solution for a flat disk. 


Explicit form of the system of equations for a twisted disk 

Now, using (12.40112. 4^ . we write equations (I2.37|) in explicit form by keeping only 
the terms in the leading order in the two small parameters 5 and u = td/feJj- 
Here, we take into account that quantities of ’thermal’ origin in the background 
solution are small, i.e. p, e, ry oc 5^p and oc oc (see Chapter 1). 

We postpone discussing the effects of the fluid non-ideality for a while. Note 
that this assumption not only corresponds to vanishing terms with viscosity 
coefficient and energy flux density, or their perturbations, but also means the 
absence of contributions oc t/'’. To select the leading-order terms in the ideal 
fluid approximation, we start by considering second terms in the r-, r- and p- 
projections of (12.371) . It turns out that such terms are proportional to 6/3 here, 
and in the r-projection of (I2.37P this contribution is due to the projection of 
the vertical pressure gradient onto the orbital plane of motion of matter in the 
twisted disk (see the analysis in Section 2.1 where this quantity was denoted 
by “(Vp)r”). In addition, in the r- and ^-projections of (I2.37|) involve terms 
oc S~^uP, which should be also kept. On the other hand, the first terms in 
the T-, r- and p-projections of (12.371) give rise to terms containing Eulerian 
velocity perturbations, as well as the Eulerian rest-mass energy density 

perturbation, pi. Hence, we conclude that 

oc max{(5, /3, and pi (xmax{S,S~^u} pfd. (2.49) 

In addition, for reasons that become clear below, we temporarily keep partial 
derivatives of w* and pi with respect to time, despite their being u~^ times 
smaller than the quantities themselves. Finally, the first terms of the r- and 
p-projections of (12.371) also contain terms with the combination d^pv^, whose 
amplitudes are restricted by the order oc max{4, (3 by equation (12.481) . 

Now, using the result (12.491) . it is easy to select the leading terms entering 
in the t-, r- and (/^-projections of (|2.37l) due to the fluid non-ideality. The 

^As we discussed above, the smallness of t^/tev is necessary to ensure that the accretion 
flow outside the equatorial plane of the black hole can be considered as a ’disk’. In turn, this 
is jointly ensured by the smallness of both 5 and t^jiLT "C 1 (see Section l2.1.1lt . 
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most troublesome here is the contribution due to the shear tensor perturbation, 
2r]af^, which appears in (see (|2.38l) and (I2.39l) i. However, most of the 
terms from this contribution contain jointly 77 oc (5^ and u® oc 6/3; therefore, it is 
clear that it is necessary to include only the terms in which the derivative with 
respect to ^ (lowering the order in S) occurs twice. This fact strongly reduces 
the number of ’viscous’ terms to be kept. Besides, by similar considerations, 
the final expressions will not contain terms with q, qi and ryi. Finally, we stress 
once again that in addition to the purely ’viscous’ terms mentioned above, the 
contribution due to the radial advection that appears in the background solution 
with non-zero viscosity should not be forgotten. We are concerned with the 
terms that can appear in the ’non-viscous’ part of the stress-energy tensor (see 
the first term in (11.591) 1 due to the non-zero value of oc <5^. 

Taking all the above into account and using relations derived in three pre¬ 
ceding Sections, we obtain the t-, r- and (/ 3 -projections of (12.371) in the form 




K3 {U^ f + {U^f 




W 


) + -U^U^d^Pi + 

J Ki r 


K^^{U^fp^+[2U^- 

+ dripU^V^) + B^pU^ + FJ = 

rd.piU^fK^U + -pU'^U^d^W, (2.50) 

Ki r 


Ko 

- 

Ki r 


Ki 


)_ I _ 1 


K, 


K 1 K 2 \ K: 


-f: = 


ZU^U^, (2.51) 


K2^rrrr^. . f ^ iU^f 


K^U^U^pi + 




TJV 

- 2ar^U‘^ ) ^pv'^F '''' ’ d^pi+2 - pd^v'^+ 

K 2 J Ki r r 


ir^ Ki K^y 

dripU^F) + d^pU^V^ + ^ ' Jj U^pv^ 

r^AiA2 


K^rd^pUW^U + ^p{U^fd^W, (2.52) 
Ki r 


where 


K = 1 — ar 




7^2 C/V ’ 

and is the total contribution due to non-zero viscous forces and the radial 

advection of matter in the background solution cx 
Explicitly, 


jjv 


F: = — {d^Tf - rWd^Tl^) - rd^pU^U'^W, F^ = 


FIf = - rWd^T^'^) - rd^pU^U''W, (2.53) 
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where 


t ;« = + U^d^W), (d^v‘^ - 2a^U^{U^fz] , 

^2 ^2 \ ^2 / 

"'"'-’"■O'- 

We note that and have the meaning of perturbations of the viscous 
stress tensor. In these expressions, the terms oc (3 contributing to the shear 
tensor perturbations appear due to the twisted basis. Conversely, relates to 
the background. Nevertheless, for the sake of brevity, we use the same notation 
with the index ‘V” for these three quantities. 

Finally, we assume in (12.50112.5^ that in the relativistic coefficients Ki, K 2 
and K 3 , the argument Rj is replaced by r, because R} = and accounting 

for the dependence on ^ here always gives rise to a small correction oc only. 

It remains to write the explicit form of the ^-projection of (I2.37p . Similarly, 
we start with the contribution of terms in the ideal fluid approximation, and 
first rearrange the first term in (12.371) . The leading-order terms in <5 here are, in 
particular, pv‘^ and pi, but additionally multiplied by This means that their 
amplitudes are restricted by the order max{d^, it}/3. Moreover, now enters 
the term “U‘^pd^v^” which also implies the raising of the order of smallness by 
S compared to (12.50112.521) (it can be seen that in formulas (12.5011 and (12.5211 
entered in combination with d^p). Besides, of all terms of a ‘thermal’ origin 
we must now keep the term with d^pi, since it also is of the order of due to 
Pi ^ S'^pi oc pS^f3. 

Turning now to the second term in the ^-projection of (I2.37L we write all 
terms through the order oc max{d^, u} /3. From similar considerations, the terms 
due to fluid non-ideality (including ‘advective’ terms proportional to oc t/'’) are 
also kept here, with their smallness increased by the coefficient S compared to 
what we did in (12.50112.521) . 

We thus obtain the following equation: 


rr</>o £ , 

+ r -h 4- 


K 2 


I — 2 ar 


K 3 \ Pi 


K2U^J p 


-+ 


'k’i 

K'2 

ar 

(KiV 


u-\] 

Ki 

K2 

2 K1K2 

Uaj 


^ ) 


K, 


dtnU — 2 a 


K 2 


= 


K:iZ , KlZ 


r K1K2 


El 

K 3 


rU^U^ + ar^§-{U^fd^U, 

Ki 


(2.55) 


where 


d^wir:^+ t:i) + {^t :^+, ( 2 . 56 ) 

and = pU'^U''. We do not provide the explicit form of here, because 
it is not required in the final form of the twist equations. 
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Everywhere in (12.55112.SBl) . except in the second term in square brackets in 
the right-hand side of (12.551) . the argument Ri in the relativistic coefficients Ki, 
K 2 and K 2 , is replaced by r. The mentioned term is an exception because this 
term alone has the zeroth order in small parameters 6 and u in equation (12.551) . 
But because we have kept the terms oc max{(5^,u} in (I2.55L in the term under 
discussion it would be necessary to take corrections oc into account due to 
the dependence of the relativistic coefficients K 2 and if 3 on We did not do 
that for the reason discussed in the next paragraph. 

2.3.2 Completing the derivation of twist equations 

Thus, we have written the twist equations in the leading orders in small param¬ 
eters 6 and u. All corrections linear in the Kerr parameter a have been taken 
into account. If we temporarily set a = 0 and consider equation (I2.55L we see 
that, on the one hand, it contains the terms proportional to the rate of change 
of the disk twist, U, and on the other hand, it has terms containing perturba¬ 
tions of the physical quantities of the order oc S^. Thus, we can say that due 
to the internal forces only, does a thin twisted disk evolve on a long timescale 
such that u ^ 6^. Then it becomes totally clear that equations (12.50112.5^ are 
restricted by the order oc S, and equation (12.551) is restricted by the order oc S^. 

At the same time, when the parameter a is non-zero, a ‘large’ term of the 
zeroth order in 6 and oc aZ arises in the right-hand side of equation (12.551) . This 
term describes the gravitomagnetic interaction of the rotating black hole with 
the tilted/twisted disk. In order that all terms in (12.551) be balanced with each 
other, we must assume that a ^ 6'^. But it then becomes clear that all additional 
corrections ^ a in equations (12.50112.5^ have the next order in S and can be 
omitted. The same relates to all terms oc ad'^ in equation (|2.55|) . including the 
correction cx 6 ^ due to the dependence of the relativistic coefficients on Rj in 
the gravitomagnetic term itself. 

In fact, this means that when considering the dynamics of a twisted thin 
accretion disk near a rotating black hole, it suffices to use the background model, 
i.e. the corresponding flat disk, in the Schwarzschild metric with a = 0. The 
assumption of the slow black hole rotation itself was needed because otherwise 
the accretion flow (including non-stationary one) could not be regarded as a 
disk, since the vertical hydrostatic equilibrium there would be violated (see 
Section 2.1). Of course, these conclusions relate to only slightly tilted/twisted 
and geometrically thin disks with /3 <C 1, <5 <C 1. 

In what follows, we therefore set a = 0 in all terms except the gravitomag¬ 
netic one. This significantly simplifies further calculations that are required for 
obtaining the twist equations in the final form. Let us first analyze equations 
(12.501) and (12.521) . It is convenient to consider their combinations which contain 
either p ov v'^. 

Eliminating v'^ for a = 0 we obtain the equation 

iU^pd^W + - rWd^r:^), (2.57) 

where we have omitted the term pi, which is of the next order in <5 compared 
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to the other terms. In the Newtonian limit as r —>■ oo, equation (I2.57P reduces 
to the continuity equation for perturbations. 

Next, eliminating pi for a = cH, we obtain the equation 


Ki 


v'^ + 


1 , 
r U'^ * 


drU'^ K[ \ 


^= o. 


(2.58) 

In the Newtonian limit, (12.581) reduces to the azimuthal component of the 
Navier-Stokes equation for perturbations. 

Finally, (I2.5ip with a = 0 takes the form 


Ki 


jjv 

- 2 




+ = IFr 


diP 


(2.59) 


In the Newtonian limit, (12.5911 reduces to the radial component of the Navier- 
Stokes equation for perturbations. 

It is important to explain why we have kept terms with and v‘^ in equations 
(12.581) and (12.591) although they are of the next order in (5. As mentioned in 
Section 2.1, in the Newtonian limit the epicyclic frequency becomes equal to 
the Keplerian circular frequency, which results in a resonance growth of the 
amplitude of velocity perturbations of gas elements in the disk under the action 
of the radial projection of the vertical pressure gradient, (Vp)^, which is limited 
only by the viscosity. Mathematically expressed, in the limit of an inviscid 
Keplerian disk, equation (I2.58|l yields in the leading order in the parameter u 
(with the term with oc omitted) such a relation between and that the 
sum of the second and the third terms in (12.591) vanishes. But, because there is 
a term oc <5/3 in the right-hand side of (12.591) . it follows that v'' (and hence v'-^ 
as well) acquires the first order in <5 in the considered case. Either viscosity or 
relativistic corrections eliminate the Keplerian resonance, and the amplitudes 
of i]'' and decrease again to the third order in <5. 

Now, from equation (12.551) we need to derive the so-called twist equation that 
plays the principal role in the twisted disk theory. For this, we need to explicitly 
determine the value d^p/p, which is done in the next Section. Although the 
Schwarzschild approximation is sufficient, we also take linear corrections in a 
into account. This is required below in obtaining an additional expression for 
the Lense-Thirring frequency in terms of the relativistic coefficients used in the 
twisted basis. 


Equation of the vertical hydrostatic equilibrium 

Let us write the ^-projection of (I2.36P in the basis Bq to the leading order in 6, 
as we did in Chapter 1 employing basis (11.49111.521) (see equation (11.741) 1. Taking 
into account that the four-velocity of the flow is {W^, 0, 0} in the leading 

order in <5, we obtain the following equation 


4 £ = 

p r 


'k'. 


K2 



K 3 K 4 

Ko 


(2.60) 


where and satisfy the normalization condition (I2.40p and the geodesic 
equation (12.441) . With this in mind, we arrive at the final form of the hydrostatic 


= 0 also in the expression for T/f^. 
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equilibrium equation 


d^P 

P 




l-2ar^^ 
K2 Uv 


(2.61) 


where the Schwarzschild profiles of t/'^ and are used in the term with the 
parameter a. 

It can be checked that with the substitution ^ —>• zjK^^ equation (12.6111 is 
equivalent to (11.7411 in the linear approximation in a. Here, we should only take 
into account that rg = ^ 2 ?’, where rs is the Schwarzschild coordinate equivalent 
to the coordinate r in (11.7411 . 


Twist equation 

Our goal is to rewrite (I2.55|l in divergent form. Without accounting for the 
gravitomagnetic term, equation (I2.55L in which we also set a = 0, must respect 
the conservation law of the angular momentum projection of the twisted disk 
onto the equatorial plane of the black hole (the conservation of the disk angular 
momentum projection onto the black hole spin in our problem, linear in /3, 
follows from equations for the background, since the corrections due to the small 
tilt are proportional to oc 1 — cos/3 ^ /3^), which reflects spherical symmetry of 
the Schwarzschild metric. 

It turns out that to do this it is necessary to eliminate from (12.551) and 
pi, in the left-hand side of this equation. Therefore, we will use equations (12.571 
I2.59|) with = i)^ = 0 for our purposes, because nowhere we will deal with 
resonance combinations of u’’ and v‘^ that vanishes in the main order in u in the 
Keplerian inviscid limit (see the comment to equations (12.5811 and (I2.59|) above). 

First, in the right-hand side of (12.5711 we rewrite the term with d^pW through 

and v'^ using (12.591) and (12.611) with a = 0. In the resulting expression for 
Pi we replace using (12.581) . Here, the derivative with respect to (p can be 
eliminated using the harmonic dependence on (p (see (12.281) ). In other words, 
dtfiip = —1- Substituting the obtained expressions for pi and in (I2.55L 
integrating over ^ and performing integration by parts wherever necessary, using 
the fact that the corresponding surface terms vanish as p —>■ 0, we arrive at the 
compact equation 



- J { dri^rKiKlU^pd^v^ + r^KiKlT^^), (2.62) 

where, as usual, E = f pd^ is the surface density of the disk, and the bar over 
means that it is integrated over In the Appendix of [5], it is shown that 
(12.621) can be used to obtain the angular momentum conservation for the twisted 
disk. 

Equations (I2.58L (12.591) and (12.621) represent a closed system of equations 
describing the dynamics of twisted configurations provided that the correspond¬ 
ing model background is specified. Unknown variables in this system include 
the velocity perturbations v'^ and v'^ and the quantity Z characterizing the disk 
geometry. We emphasize that in deriving these equations we essentially used 
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only three main assumptions: a <C 1, <5 1 and /3 <C 1. This means that 

the equations describe the dynamics of any geometrically thin accretion flow 
(disk) with any parametrization of viscosity, any radial and vertical structure in 
both stationary and non-stationary case. In the latter case, we mean the non¬ 
stationary background: the equations determine not only the dynamics of twist 
perturbations propagating in a stationary flat disk, but also the dynamics of 
twisted rings/tori, when the evolution of the geometrical form occurs in parallel 
with its expansion in the radial direction due to turbulent viscosity, which also 
results in the evolution of the background itself. 


Once again about the characteristic frequencies of the problem 

In Section 12. 1. 11 we already obtained relativistic expressions for the charac¬ 
teristic frequencies of the problem. These include the circular and epicyclic 
frequencies of free equatorial motion, as well as the frequency of vertical oscil¬ 
lations and the precession frequency of tilted orbits. Here, we wish to obtain 
expressions for these frequencies, but now in terms of the values used above to 
construct the theory of twisted disks, i.e. in the basis Bq. These expressions 
are required to write the twist equations in a more compact form. 

The circular frequency of the free equatorial motion as measured by the clock 
of an infinitely remote observer, which we already presented in equation (I1.16L 
can be obtained simply by dividing d(p/ds by dr/ds given in (12.471) . We obtain 


n = 


Ki 

KK2^' 


(2.63) 


Using (|2.45l) and (12.401) . and also remembering that rg = rK 2 , we can check 
that (I2.63|l coincides with (I1.16|l in the linear approximation in a. 

We now consider small vertical deviations from the circular equatorial mo¬ 
tion. In Section !^. 1. 11 we discussed that the frequency of vertical oscillations as 
measured by an infinitely remote observer, is the locally measured frequency, 
divided by the t-component of the four-velocity of circular motion, 11*. The 
frequency explicitly enters the equation of hydrostatic equilibrium (see equa¬ 
tion (ll.74|) or equivalent equation (12.6111 with the substitution ^ 0 / 1 ^ 2 )■ Using 

relations (I2.47L we express U* = dr/ds in terms of U'^: 


Ul = KK^^U\ 


whence 




KU'^ 


n,: 


Ki rK2 VL ’ 
where the final expression was obtained using (12.6311 . 

But then, from a comparison of (12.6411 with (I2.61L we see that 


(2.64) 


— H 


. ^3 u- 

1 — ar - 

K2 Uv 


(2.65) 


where the Schwarzschild profiles for and are used in the term with the 
parameter a. 

Then, using (12.1211 . we obtain the Lense-Thirring frequency 

nLT = a^^. ( 2 . 66 ) 
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It suffices for our purposes to know the epicyclic frequency in the Schwarzschild 
case with a = 0. This expression can be the most easily derived directly from 
the twist equations, more precisely, from their part that describes the dynamics 
in the plane of disk rings, i.e. from (12.5811 and (12.591) . Setting the ‘viscous’ terms 
and radial projection of the pressure gradient in the right-hand side of (12.5911 
equal to zero, as well as omitting the dependence of v'" and on </?, we obtain 
equations for the Eulerian perturbations that describe a free motion of gas ele¬ 
ments slightly deviating from the circular motion. Clearly, these equations are 
equivalent to (I2.5ll2.(i|l which were written in the basis (11.49111.5^ . From these 
equations, we obtain the following equation for u’': 


+ 2 


KiK{ 


fdrU^ Ki iun^ \ 
\ U'fi ^ Ki [U^Y J 


= 0 , 


(2.67) 


where the expression before u’’ is equal to It can be rewritten in a more 
compact form 

TJT\I 

( 2 . 68 ) 


^2^^K[{K,Un' 


to ensure that it coincides with (12.9L considering that the radial Schwarzschild 
coordinate rs = rK 2 enters the last equation. 

Finally, for convenience, we introduce the following quantity with the di¬ 
mension of frequency that appears in our problem. In the Schwarzschild case, 
a = 0, 

K[ 1 rs-3 


n = 


K2U^Uv r|(rs-2)1/2’ 


(2.69) 


which tends to the Keplerian value in the Newtonian limit. 

Using (12.6811 . (12.6811 and (12.6911 allows us to write equations (I2.58|l and ()2.59ll 
in a more compact form. Lense-Thirring frequency (I2.66|l . evidently, enters the 
gravitomagnetic term in (12.6211 . However, we deal with this rewriting in the 
next Section when considering a specific background model. 


2.3.3 Twisted equations in the particular case of a sta¬ 
tionary vertically isothermal a-disk 

We now consider the form the twist equations take in the specific background of 
a stationary a-disk which we have discussed in Chapter 1. This does not mean, 
however, that only stationary twisted solutions are to be considered. In other 
words, the equations we obtain are also applicable to arbitrary non-stationary 
dynamics of the corresponding twisted disk. For example, they enable us to 
calculate the evolution of the shape of an (infinite) initially flat disk instantly 
tilted to the equatorial plane of a rotating black hole. The initial stage of the 
evolution of such a disk was qualitatively described in Section 12.11 In addition, 
these equations describe the wave-like (in the case of a disk with sufficiently 
small a < 5] see also [33]) or diffusion-like (in the case of a disk with sufficiently 
large a > 5\ see also [33]) dynamics of some twisted perturbation imposed on 
the disk lying initially in the equatorial plane of the black hole. 

Explicit form of the necessary background profiles 

The twist equations contain the quantity (as well as ry), related to the corre¬ 
sponding flat disk model. We could obtain the explicit form of these quantities 
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by integrating the t— and (^—projections of equations (12.361) . However, it is 
simpler to use the results of Chapter 1, where we have already obtained this 
quantity, denoted by there (see equation (11.901) ). We should only take into 
account that now we are working in another basis than that used for the flat 
disk, and therefore the transition from to should be specified. First, 
using the orthogonality condition for the shear tensor and, hence, for the vis¬ 
cous stress tensor, (ll.66|) . we see that only one component of the viscous stress 
tensor, T^,^', is non-zero in basis (jl.49m.5^ . because the four-velocity there has 
only non-zero time component up to the terms oc The prime here marks 
basis (|1.49111.5^ . Further, the (orthonormal) bases are different only in that 
the observers associated with basis (11.49111.5^ move in the azimuthal direction 
with the velocity of the free equatorial circular motion, whereas the basis Bq 
corresponds to the observers at rest. Therefore, the transformations of vectors 
and tensors must be equivalent to the usual Lorentz transformations. Using |35] 
(see exercise 1, paragraph 6 therein) we see that where U'^ is 

the Lorentz factor of azimuthal motion. Finally, it should be additionally taken 
into account that integration over ^ differs from that over z by the coefficient 
K 2 As a result, we obtain 

T7 = (2.70) 

A 2 

We note that it is possible to change from T^'^' to using the relation (11.341) by 
writing it for two bases, equating the right-hand sides and then multiplying one 
of the sides of the obtained equalities by matrices inverse to the basis matrices 
there. Here, we should only take into account that in the basis Bq the radial 
coordinate was changed, (12.141) . i.e. that rs = rK 2 in the notation of this part 
of the paper. 

Next, for the case a = 0, which is sufficient here, it is easy to express Ti, in 
terms of elementary functions. Indeed, the integral in (11.901) can be taken by 
the substitution y = ^Jrs'- 


E 


1/2 


j f y 1 

C J y^-5 2 


For T'^‘^ with account for (I2.70|) we finally obtain 


y + Vs 

y-V^ 


M 

T^v = —u^r- 
27r 


•3/2. 


L{r) 


kTkI' 


(2.71) 


where 


i = 1 _ ^ ^ In " y3)(3 + 2^) 


2 y 


{y + 'Ji) 


(2.72) 


As it must be, L = 0 at rs = 6. We note that L = Y{a = Q) where Y was 
defined in ()1.101|) . 

On the other hand, the expression for in (11.721) . in our case a = 0 in the 
basis Bq, can be rewritten in the form 


tr ^ - 


D 


3 / 2 ^ 


[/' = ^Klu^iru- = ~unu-) 


where, as usual, we use relations (12.471) . 
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Then 


(2.73) 


f^V = If, - LlTjrffJV_ 

2 


As in Section fl.3.61 equating expressions (12.7111 and (12.731) we obtain 


M 


r.-l/2 


L 


V = 


37r ' 


(2.74) 


In the Newtonian limit, far away from the inner edge of the disk, equation (12.741) 
gives a well-known result fj = M/ (37r). 

We assume that the kinematic viscosity is proportional to the characteristic 
disk half-thickness times the sound velocity in the disk: 


cyCshrj 


(2.75) 


where hp is the proper characteristic disk half-thickness, which in our coordinate 
system is hproper = K 2 h and a is the Shakura parameter, which is assumed to 
be constant. Because (12.611) implies that Cs ^ \/P/P ^ U'^^h/r, we finally define 
a by the equality 

= aK2U^h?/r. (2.76) 


Using (12.741) and (I2.76L we obtain the relation 


= 


M 


r.i/2 


, —V 

37ra \^U'^(17‘^)2 ^3^5/2y 


(2.77) 


To find U’' in the advective term in (I2.62|l . we use the rest-energy conser¬ 
vation law in the basis Bq for the stationary disk. Again, we use result (11.861) . 
Recalling the transition to the isotropic radial coordinate, the relation between 
the coordinate and physical velocities (11.471) and (I2.47L and the difference in 
the definitions of E, we obtain 


Then [/’’ can be derived from (12.7811 and (12.771) as 

3a 




(2.78) 


(2.79) 


Finally, we need to know the profile 5{r). Note that this value is invariant 
under the transition between the bases ()1.49111.5^ and Bq, since the change 
from hp to h and from rg to r is scaled with the same coefficient K 2 . 

In the gas-pressure-dominated disk with the Thomson scattering opacity, it 
follows from (11.1031) with a = 0 that 




(2.80) 


In order to derive a simpler form of twist equations, we need to specify the 
vertical profile of the rest-energy density. Here we use its simplest form in an 
isothermal disk: 


p = Pc exp 

where pdr) is the equatorial density. 


2/i2 


(2.81) 
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Transition to complex amplitudes 

In the case of an isothermal disk, the velocity perturbations u’’ and v'^ taken in 
the form 


= ^(sin (/j + ^2 cos 1 ^) z)’’= ^(i?i sin </? + i ?2 cos (^) (2.82) 

satisfy equations (12.581) and (12.5911 provided that v does not change with the 
height and the amplitudes Ai, A 2 , Bi and B 2 are functions of r and r. Indeed, 
in this case, all ‘thermal’ terms are oc and the dependence on ^ with the 
ansatz (|2.82|1 is identically satisfied in the considered equations. 

Let us introduce the complex amplitudes 

A = A 2 +iAi, B = B 2 + iBi and W ^ + i ^2 = (2.83) 

By composing two combinations, (I2.58P + i&n (|2.58l) and (12.5911 + i9, ., (12.591) . 
we see that all terms in these combinations are oc e”*"^. In particular, the terms 
containing W and d^pW pass into the terms respectively containing —iW'e”*''” 
and AV'e"*"^. 

As a result, we obtain the following complex equations 

^2 o 

A-[i- a)nA + = --iaKiiW^fU^nW, (2.84) 

B-{i- a)nB - 2nA = -{i + a)U^Vt^A^', (2.85) 

where we have used equation (|2.76|) . as well as expressions for frequencies (12.6311 . 
(12.681) and (|2.69ll obtained in Section [2.3.21 

In a similar way, by using (12.821) and (12.831) and composing the combination 
(12.621) +ic},., (l2.62D . we derive an equation for complex amplitudes. In the right- 
hand side of this equation, the integration over ^ should be performed under the 
derivative with respect to r. For an isothermal disk with density distribution 
(I2.8II1 . the equality / holds. Thus, the derivative with respect 

to r acts on terms proportional to or fj. Instead of these combinations, 
we substitute equations (12.771) and (12.741) there and group common constant 
factors before the derivative with respect to r. Additionally, instead of and 
we substitute expressions (|2.79|1 and (12.7111 in the left-hand side of the 
discussed equation and then divide the whole equation by S. The obtained 
equation contains M and E only in the combination M/T,, which we express 
through (5^ and other known quantities using (12.7711 . Also using the expression 
for Lense-Thirring frequency (12.661) . we finally arrive at the following equation 


3 K'^ 

W - tfliTW -t - Ki{rK2) 

2 K 2 




S^KfU^ d 
2A/^kI/^L dr 


.3I2k^/2 




jjv 


( {i + a)B + aU^W ) 


( 2 . 86 ) 


Equations (12.84112.SHIl form a closed system of equations for the quantities 
A, B and W as functions of r and r. In the weak gravity limit they reduce to 
equations (30), (31) and (33) in [36] . 


50 



















































2.4 Stationary twisted disk 


2.4.1 Main equation and boundary condition 

We now consider stationary solutions of the system of equations 
The main goal of this Section is to calculate the shape of a stationary twisted 
disk. 

We set A = B = W = 0. After eliminating A from (12.84112.85]) we obtain 


(* - a)2q2 


1 + 7- {i + a)U^VL - KAU^fU^n W', 

L I — a \ 

(2.87) 

whence we express B through W' and substitute it in equation (12.861) . We thus 
obtain the equation 


Ki 


_(liT 

\KiU^ 


dW 


i,dW 


r{a,rs)— -3aU^{l-L-^) — 


4ia 


drs 


drs 6'^KfrgU‘^ 


where the asterisk denotes complex conjugation and 


/(a, rs) = (1 + - 3iaKf) 


rs(» - a) 

ars{a + 2i) — 6 


+ a. 


W = 0, 

( 2 . 88 ) 

(2.89) 


We note that equation (12.881) was written after changing to the Schwarzschild 
radial coordinate rg. In what follows, we wish to consider only the case a > 0, 
i.e. a prograde disk. It can be seen that the problem has two free parameters. 
First of all, this is the combination 6 = i5*/-\/[a[. Clearly, 6 ranges from 0 to 
oo and characterizes the relative role of the hydrodynamic and gravitomagnetic 
forces acting on the disk rings. Second, (12.881) contains the disk viscosity param¬ 
eter 0 < a < 1. Equation (I2.88|) in the rigorous Newtonian limit with non-zero 
viscosity reproduces the corresponding equation (2.10) from [37] and, addition¬ 
ally, with post-Newtonian corrections, reproduces equation (33) from [25j . what 
is checked in [^ (see paragraph 4.1 therein). 

The coefficients of equation (12.881) have a singular point at the inner edge 
of the disk at rs = rs = 6, where L vanishes. The regularity of the solution 
at rs must yield a condition for the function W. Using this condition as the 
initial one, we can integrate (12.881) from rs to infinity and to obtain the form 
of the stationary twisted disk. We expand equation (12.881) in series in orders of 
the small Xq = rs — rs ^ In practice, to do this, all quantities that take 
non-zero values at rs should be set exactly equal to these values and L should 
be expanded to the main order in xq. From (12.721) . we hnd 


(2.90) 

whence we see that another quantity in (12.881) that vanishes at the inner disk 
edge, (5, can be written as 

d — dmsXQ , 

where e is the power-law exponent L in equation (12.801) . Accordingly, Sms is also 
given by equation (I2.80|) , which is taken at rs and into which we now substitute 
the coefficient 72“^ from ()2.90l) instead of L. 
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After that, it is easy to obtain the equation valid for xq 1, 


dxo \ ° dxo 


Cixl-^^W + C2^=0, 

dxQ 


(2.91) 


where 


and 


2 i flLT 

f{a,rs) Uv KfrsS^, 

216a 

f{a,rs) rs 


are taken at rs. We see that for any finite viscosity the last term in (12.911) 
becomes dominant sufficiently close to the disk edge; therefore, the boundary 
condition can be straightforwardly written as 


dW 

dxo 


= 0 . 
rs 


(2.92) 


On the other hand, from (I2.9ip with a = 0, we obtain a simpler equation 
whose solution is a Bessel function: 


W = Cx~^^‘^J^^{z), (2.93) 

where 

As ccq —> 0 (12.931) tends to a non-zero constant but with a zero derivative with 
respect to xq. Therefore, in this case we return to condition (12.921) . 

Due to the linearity of the problem, it suffices to take an arbitrary non-zero 
value of W in fs, to set the first derivative of W at fs equal to zero, and with 
these boundary conditions to integrate (12.881) up to infinity. The modulus and 
phase of W give the profiles f3{rs) and 7 (^ 5 ) for a stationary twisted disk. In 
what follows, we normalize the profile /3 to unity at infinity. 


2.4.2 Disk with a marginally small viscosity 

We consider a disk with a very low viscosity separately. Clearly, it is possible to 
analytically treat the accretion disk setting formally a —>■ 0 , if simultaneously 
M —>■ 0. In such a disk, U'' —>■ 0, however, it then has definite profiles of E and 
h. 

In addition, to obtain an analytical solution, we consider the case <5 <C 1; in 
other words, we assume a sufficiently thin disk around a rapidly rotating black 
hole. 

Setting a = 0 in (12.881) yields 


+AW = 0, 

drs V drs ) 

where 

^ 24aL 


(2.95) 


(2.96) 
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The coefficients in (12.951) take real values; therefore, there are real solutions 
of this equation. This means that in the absence of viscosity in the stationary 
twisted disk 7 = const which can be set equal to zero by the corresponding 
choice of the reference frame. Therefore, the variable W is identical to the 
angle /3 in this Section. 


The form of the disk near its inner edge 


In the foregoing, we have already presented the solution near the inner edge of 
the inviscid disk (see equation (12.931) ). The constant Ci in this case has the 
explicit form: 


Cl 


24aC/^ 


(2.97) 


and is taken at rs = rs- 

Using the known approximation for the Bessel function of small argument, 
we obtain a relation between the constant C in (|2.93|) and the value of W at 
rs, W(rs) = Wo: 


C' = r 


/ 3-4e 

V2(l-2e) 


Vx \ 

2{l-2e)) 


-l/2(l-2e) 

Wo, 


(2.98) 


where r(x) is the gamma-function. 

In addition, we need the asymptotic of (12.931) for z ^ I. Clearly, z can be 
large, even for a: <C 1, because ^ 1. Hence, for z » 1 we obtain 


W 



TT 1 — e \ 

2 1 - 2 ej ■ 


(2.99) 


The form of the disk at large distance 

We consider equation (12.881) for rs » 1 and a —>■ 0. Importantly, we can¬ 
not set all variables in (I2.88P to their Newtonian values and have the vis¬ 
cosity simultaneously vanished. This already follows from the fact that then 
f{a,R) —>■ 1 /( 20 ) —>■ 00 . Physically, this reflects the fact that, as we mentioned 
above, in the absence of viscosity in the strictly Newtonian potential, a Keple- 
rian resonance occurs when the circular and epicyclic frequencies coincide, and 
perturbations in the twisted disk grow infinitely due to the action of the radial 
projection of the vertical pressure gradient. Therefore, a stationary twist is im¬ 
possible in this case. Taking the next-order term in the expansion of f{a,rs) 
in small into account, we obtain 


f{oi,rs) 



( 2 . 100 ) 


As a —> 0, f{oi,rs) now remains finite at any finite rg. Nevertheless, it makes 

the leading contribution due to relativistic effects, and all other variables in 

_ 1/2 

(I2.88|) can now be set equal to their Newtonian values = 1, 17'^ = j 
L = 1 and Ki = 1. Moreover, we neglect the weak dependence of 5 on rs far 
from the black hole and set 6 = 6^. 
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_ 1 /o 

After that, by introducing the new independent variable xi = Vg ' -Cl, 
we obtain the equation 


^2 j 

- 2—W + 96(5"2a;fW = 0. 
dx^ dxi 


( 2 . 101 ) 

The solution of (12.10111 can again be expressed in terms of a Bessel function: 

W = J_3/5(zi) + A2 J3/5(^l)), (2.102) 


where 


Zl 




(2.103) 


Ai and A 2 are constants. 

When rs is so large that <C 1, the first and second terms in (I2.102L 
multiplied by x-{ , respectively tend to a non-zero constant and to zero. This 
allows us to express the constant A\ in terms of the value of W at infinity, Woo: 


W — 

VV on — 


/ 5 x3/5 ^3/5 

Uvej r(2/5) '• 


(2.104) 


In the opposite case, at zi » 1, i.e. closer to the black hole, we obtain 
another asymptotic form: 


W 




.,- 1/8 


27rv^ 


Ai cos -f A 2 sin . (2.105) 


WKB-solution for the shape of the disk 

Throughout the disk the asymptotic solutions (j2.99ll and (12.10511 can be matched 
by a WKB-solution of equation (12.9511 . Indeed, because we are considering the 
case 5 <C 1, the ratio of A and b in ()2.95|) . A = A/ 6 , is large at all rs such that 
z and zi are large. 

The WKB-solution has the form 

W ft! cos y^drs + (pwKBj^ , (2.106) 

where the constants C 3 and (j)WKBj should be chosen such that (12.10611 is 
smoothly matched with formula (12.9911 in the corresponding region. It can be 
checked that this yields 

(j^WKBj = i_2e (2.107) 

and 

where we assume that Ki and are taken at rg = fg = 6 and L k, 112 near 

rs- 

Next, in the limit rs ^ 00 we can set A and 6 before the cosine in (12.10611 
equal to their Newtonian values. In addition, the integral in (12.10611 can be 

represented as I{rs) = \/xdrs = I — ^/xdrs, where / = \/\drs- 
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Taking into account that the Newtonian value A = 24<5 we have that 

Vxdrs ~ and therefore 

W « ^ 3 ^cos -I- cPwKB^j ■ (2.109) 

Solution (12.1091) must be smoothly matched with expression (12.1051) in the cor¬ 
responding region, which yields the values of constants Ai and A 2 . It can be 
checked that they are 

= )J~^C3 cos (^I + (j)wKBj - j *^ 0 ® Yo ’ 

A 2 = y^Casin + (I>wkbj + ^)/ (2.110) 

Thus, equations (12.931) . (12.1061) and (12.1021) jointly with coefficients (12.1081) . 
(12.1101) and phase (12.1071) . determine the shape of an inviscid stationary rela¬ 
tivistic twisted disk at all distances in the range from rg = fs to rg = 00 . 


Resonance solutions in the low-viscosity disk 


We note that equations (I2.98L (I2.108L (12.1101) and (12.1041) provide a relation 
between Wq and Woo: 


Woo = aot(5)Wo, 


( 2 . 111 ) 


where the explicit form of Ctot{5) follows from these formulas. In particular, as 
follows from (12.1041) and (I2.110L Ctot{S) oc cos(/ -I- 4’wkbj ~ ^)- 

We hence conclude that for some discrete set of S for which cos{I + (j)w kb j — 
^) = 0 so Woo = 0 despite that Wq ^ 0 . 

From equations (12.961) . it is possible to represent the integral I in the form 
/ = S~^I, where I does not depend on S. This allows us to write the singular 
values of S explicitly: 


4 


/ 

§(w + ^+2k)’ 


( 2 . 112 ) 


where k is an integer number. 

The values of Sk correspond to such a balance between the external grav- 
itomagnetic force and the internal pressure gradient in the disk that leads to 
the disk twist even if the matter flowing into the disk at infinity moves in the 
equatorial plane of the black hole. We note that the solution in the form of a flat 
disk lying entirely in the black hole equatorial plane, of course, also exists for 
these 4- This non-uniqueness of the solution disappears for any small viscosity 
in the disk, for which Woo = 0 always implies Wq = 0. For small a <C 1, the 
disk ’feels’ these ’resonance’ solutions, and its inner parts deviate significantly 
from the equatorial plane of the black hole, even when the outer parts of the disk 
lie almost in the equatorial plane. Figure shows the curve corresponding 
to the analytical solution ( 12 . 1111 ) , as well as several curves for a viscous twisted 
disk obtained by integrating the original equation (12.881) . We see that already 
for a = 10 “^ the discussed resonances are almost entirely suppressed. 
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Figure 2.1: Ratio of the tilt angle of the inner disk edge to the tilt at infinity, 
/3o//3oo) as a function of the parameter S. The solid curve shows the numerical 
solution of equation (12.881) with a = 0, the dotted curve represents the analytical 
dependence where Ctot is given by equation (12.1111) . The dashed, dash- 

dotted and dash-dot-dashed curves are obtained by numerical integration of 
equation (12.881) with a = 10“"^, 10“^ and 10“^, respectively. 
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Figure 2.2: Contours of constant ratio Po/f3oo on the parameter plane {S, a). 
The numbers show the value of Po/Poo for each curve. The dashed curve in 
the right part of the Figure separates the region where the change in /3 with rg 
is more than 10% of /3oo (to the left) from the region where the disk twist is 
insignificant, and /3 deviates by less than 10% from jSoc (to the right). 

2.4.3 The disk behavior on the plane of parameters a and 

(5 

In conclusion, we present the full study of regimes of behavior of a stationary 
twisted relativistic disc near a rotating black hole. It is convenient to show 
the results of numerical integration of equation (12.8811 on the plane of free pa¬ 
rameters of the problem, 5 and a. The first parameter varies in the range 
10“^ < 5 < 10 and the second parameter varies in the range 0 < a < 1. As 
follows from Fig. (lO) . at small 5, i.e. when the gravitomagnetic force exceeds 
the internal forces in a twisted disk, it either lies in the equatorial plane of the 
black hole, i.e. /3o//3oo —>■ 0, or, conversely, the tilt of its rings strongly increases 
in the inner parts of the disk, with oscillations of I3{rs) along the radial coordi¬ 
nate. Note that for small viscosity, these oscillations become so strong that the 
corresponding gradient of the tilt angle, (3', leads to supersonic perturbations of 
the velocity components, Vr and v'^, at heights of the order of the disk thickness, 
^ h. This, in turn, must lead to the generation of various hydrodynamic in¬ 
stabilities and sound waves, which cause additional disk heating (and hence also 
an increase in i5), as well as the growth of a. These processes should partially 
suppress the oscillations of f3 discussed above. 

The disk alignment into the equatorial plane of the black hole occurs at 
sufficiently high viscosity, when the condition a > (5 is satisfied with a large 
margin, and is referred to as the Bardeen-Petterson effect [38]. It is seen from 
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Figure 2.3: The dependence of /3 on rs along the curve in Fig. (23) for which 
Po/Poo = 1. 5 takes the values 10“^, 10“^, 10”^ and 1 for the solid, dashed, 
dash-dotted and dash-dot-dashed curves, respectively. 


Fig. (12.211 that this effect occurs only in sufficiently viscous and thin disks. 
But already for 5 ^ a the ratio Po/Poc becomes of the order of unity, which 
means the absence of the disk alignment. At the same time, oscillations of P 
disappear. Figure (12.31) shows the profiles of P{rs) when PojPoo = 1 for several 
5. It is seen that for not very small 6, the twisted disk has a sufficiently smooth 
shape, which suggests the possibility of the existence of such configurations in 
nature. We note that P behaves non-monotonically: it first decreases and then 
increases with the decrease of rg. The latter can have important implications 
both for the disk structure itself and and for its observational manifestations. 
For example, the hot inner regions of such a disk should illuminate its outer 
parts much stronger compared to the flat disk case. Clearly, this is due to the 
disk inner parts being tilted with respect to the outer parts. 

In the region where S is of the order of or greater than unity, the action of 
the gravitomagnetic force becomes insignificant, and the disk is weakly twisted. 
In Fig. (1^ . the area to the right of the dashed line is where P{rs) deviates 
from Poo by less than 10%. It is also worth noting that for 5 > O.I the Bardeen- 
Petterson effect is completely absent for any a. 

2.5 Conclusion 

We have presented a detailed technical derivation of the governing equations for 
the evolution of the shape of a relativistic twisted disk, as well as for perturba¬ 
tions of the velocity and density inside it. Only three simplifying assumptions 
have been used: the smallness of the disk aspect ratio, ^ <C 1, the slowness 
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of the black hole rotation, a <C 1, and the smallness of the disk’s rings tilt to 
the equatorial plane of the black hole, ,5 <C 1. This allowed us to formulate 
equations (|2.58l) . (I2.591) and (12.621) for three variables describing Eulerian per¬ 
turbations of the azimuthal velocity, v"^ and and the geometrical form of the 
disk, Z. In general, the dependence of v'^ and v‘^ on the twisted coordinates r, 
^ and T, and the dependence of ^ on r and r should be found. In accordance 
with equation (12.281) . all these variables depend harmonically on the azimuthal 
coordinate. The governing equations contain the profiles of the background so¬ 
lution, representing an accretion disk with similar radial and vertical structure 
but lying in the equatorial plane of the black hole. We note once again that 
not only the twisted disk but also the background itself can be non-stationary, 
because when deriving the system of equations (12.58112.59112.621) , the only one 
assumption about the background, the smallness of i5 <C 1, was used. There¬ 
fore, the twist equations enable us also to study the evolution of tilted/twisted 
gaseous tori/rings near the rotating black holes as they are spreading in the 
radial direction, in other words, as the non-stationary accretion goes on due to 
turbulent viscosity. 

In the particular case of a stationary, vertically isothermal background with 
a-parametrization of the viscosity, the twist equations have been reduced to 
simpler equations (I2.84|l . (12.851) and (I2.86|l for the complex amplitudes A and 
B describing the velocity perturbations, and W describing the disk geometry, 
which depend only on r and r. Here, the solution for a flat relativistic disk, 
which was presented in detail in Chapter 1, was utilized. The corresponding 
stationary problem can be described by a second-order linear differential equa¬ 
tion for W (see equation (12.881) ). The analytic integration of this equation for 
a formally inviscid disk with i5 <C 1 enabled us to find the singular resonance 
solutions for a discrete set Sk, which in fact correspond to an instability of a flat 
non-tilted disk, when the latter can acquire a twisted shape near the black hole, 
even with its outer part lying in the black hole equatorial plane. This instabil¬ 
ity, however, rapidly disappears already for a ^ 10“^ and for a > S provided 
that <5 <0.1, the numerical calculations show the Bardeen-Petterson effect. At 
the same time, already for a ^ S, the alignment of the inner parts of the disk 
into the equatorial plane of the black hole is absent, and for 5 > 0.1 smooth 
but non-monotonic profiles P{r) appear (see Fig. (12.31) ), which suggests their 
stability under perturbations and the possibility of their realization in nature. 
The last effect is confirmed by the first numerical simulations of tilted thin rel¬ 
ativistic accretion disks with J^a^a^O.l carried out in recent papers |39j 
and |40j . In these papers, a comparison with the semi-analytic model based on 
the solution of the system of equations (12.58112.59112.621) was also done for a 
slightly tilted vertically barotropic torus. 

Observational confirmations of the existence of twisted accretion disks around 
rotating black holes have just started emerging. Apparently, one of the most 
direct pieces of evidence of their existence is the observation of maser sources 
at subparsec scales in the disk around a supermassive black hole in the nucleus 
of NGC 4258 [H], [25]. The subsequent modeling in [23] and [22] showed that 
the disk twist in this case can be due to the Bardeen-Petterson effect. In the 
recent paper m, observations of jets in the nucleus of NGC 4248 were used 
to independently estimate the black hole Kerr parameter a ^ 0.7 and, in a 
similar model, to calculate the radius of the disk alignment into the equatorial 
plane of the black hole in agreement with observations. Additional but more 
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indirect arguments favoring the presence of twisted disks in galactic nuclei were 
obtained, for example, in [?5] and 113, where the observed profiles of X-ray iron 
line Ka were calculated for different accretion disk models. It was concluded 
that in many cases, the observed line profile can be more easily explained in the 
model of the twisted disk than in the model of a flat disk but, for example, with 
specific radial intensity distribution. In |48j . a similar modeling of hydrogen 
Balmer lines was performed, which should arise due to the heating of the outer 
parts of a twisted disk by hard emission from its inner parts, which have much 
larger tilt than in the case of a flat disk. The presence of twisted disks is also 
suspected in binary stellar systems with black holes. For example, this can be 
the case in two microquasars, GROJ1655-40 and V464I Sgr, in which the tilt 
of jets relative to the orbital plane was discovered (see [33] and [50]'). 

As mentioned above, equations (12.581) . (12.591) and (12.621) also describe the 
non-stationary dynamics of a torus tilted to the equatorial plane of a black 
hole. If (5 > a, the action of the gravitomagnetic force must lead to a solid-body 
precession of the torus, because in this case the twisted waves propagating at 
almost the speed of sound smear out the dependence of 7 on r due to the 
Lense-Thirring effect. Similar non-stationary models are invoked to explain the 
variability of Balmer line profiles, as well as the precession of jets in active nuclei 
(see, e.g., ED). In many papers, the precessing tori are used to explain low- 
frequency quasi-periodic oscillations in X-ray binary systems (see, e.g., [52]). Of 
special interest is the modeling of observational appearances of a tilted accretion 
disk around the black hole in the center of our Galaxy [53]. 

The theory of relativistic twisted disks presented here can also be successfully 
applied both to constructing self-consistent models of individual objects and to 
making further theoretical predictions on the dynamics of accretion flows around 
rotating black holes. 
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